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ABSTRACT 

Mass modelling of spherical systems through internal kinematics is hampered by the mass 
/ velocity anisotropy degeneracy inherent in the Jeans equation, as well as the lack of tech- 
niques that are both fast and adaptable to realistic systems. A new fast method, called MAM- 
POSSt, is developed and thoroughly tested. MAMPOSSt performs a maximum likelihood fit 
of the distribution of observed tracers in projected phase space (projected radius and line- 
of-sight velocity). As in other methods, MAMPOSSt assumes a shape for the gravitational 
potential (or equivalently the total mass profile). However, instead of postulating a shape for 
the distribution function in terms of energy and angular momentum, or supposing Gaussian 
line-of-sight velocity distributions, MAMPOSSt assumes a velocity anisotropy profile and 
a shape for the three-dimensional velocity distribution. The formalism is presented for the 
case of a Gaussian 3D velocity distribution. In contrast to most methods based on moments, 
MAMPOSSt requires no binning, differentiation, nor extrapolation of the observables. Tests 
on cluster-mass haloes from ACDM dissipationless cosmological simulations indicate that, 
with 500 tracers, MAMPOSSt is able to jointly recover the virial radius, tracer scale radius, 
dark matter scale radius and outer or constant velocity anisotropy with small bias (<10% on 
scale radii and <2% on the two other quantities) and inefficiencies of 10%, 27%, 48% and 
20%, respectively. MAMPOSSt does not perform better when some parameters are frozen, 
and even particularly worse when the virial radius is set to its true value, which appears to 
be the consequence of halo triaxiality. The accuracy of MAMPOSSt depends weakly on the 
adopted interloper removal scheme, including an efficient iterative Bayesian scheme that we 
introduce here, which can directly obtain the virial radius with as good precision as MAM- 
POSSt. Additional tests are made on the number of tracers, the stacking of haloes, the chosen 
aperture, and the density and velocity anisotropy models. Our tests show that MAMPOSSt 
with Gaussian 3D velocities is very competitive with other methods that are either currently 
restricted to constant velocity anisotropy or 3 orders of magnitude slower. These tests sug- 
gest that MAMPOSSt can be a very powerful and rapid method for the mass and anisotropy 
modeling of systems such as clusters and groups of galaxies, elliptical and dwarf spheroidal 
galaxies. 

Key words: methods: analytical - galaxies: kinematics and dynamics - galaxies: haloes - 
galaxies: clusters: general - dark matter 



1 INTRODUCTION 

The determination of mass profiles is one of the fundamental issues 
of astronomy. Subtracting the mass density profile of the visible 
component, one deduces the dark matter (hereafter, DM) density 
profile, which can be confronted to the predictions from cosmolog- 
ical iV-body simulations. This is especially relevant given the dif- 
ferences between the total NFW (Navarro, Frenk, & White 1996) or 



better Einasto (Navarro et al. 2004) density profiles derived in dis- 
sipationless simulations of a single dark matter component on one 
hand, and the 1/r 2 density profiles found for the DM in hydrody- 
namical cosmological simulations (Gnedin et al. 2004). Moreover, 
the knowledge of the total density profile serves as a fundamen- 
tal reference, relative to which one can scale various astronomi- 
cal tracers such as the mass density profiles of the stellar, gas and 
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dust components, as well as the luminosity in different wavebands. 
These studies can be performed as a function of system mass and 
other attributes such as galaxy colour (e.g., Wojtak & Mamon 2013 
and references therein). 

Mass profiles can be derived from internal motions, or alter- 
natively from X-ray or lensing observations. This paper focuses on 
mass profiles from internal kinematics. In this class of mass mod- 
eling, one has to deal with a degeneracy between the unknown ra- 
dial profiles of total mass and of the velocity anisotropy (hereafter 
'anisotropy') 



/3(r) = 1- 



(1) 



(note that, in spherical symmetry, one must have cr^ = ag). While 
radial outer orbits are expected for structures in an expanding uni- 
verse (e.g. Ascasibar & Gottlober 2008 for dark matter particles 
and Ludlow et al. 2009 for subhaloes), the dissipative nature of 
the gas dynamics is expected to produce tangential orbits in the in- 
ner regions of systems formed from gas-rich mergers or collapse. 
Therefore, lifting the Mass - Anisotropy Degeneracy can provide 
useful constraints on the formation of the structure under study. 

A common method to extract the mass profile is to assume that 
the line-of-sight (hereafter, LOS) velocity distribution, at given pro- 
jected radius, is Gaussian (Strigari et al. 2008; Battaglia et al. 2008; 
Wolf et al. 2010). These methods perform adequately on the mass 
profile, but provide weak constraints on the anisotropy (Walker 
et al. 2009). Merritt (1987) pointed out that anisotropic models have 
non-Gaussian LOS velocity distributions. Therefore, the observed 
kurtosis of the distribution of LOS velocities serves as a powerful 
constraint to the anisotropy (Gerhard 1993; van der Marel & Franx 
1993; Zabludoff, Franx, & Geller 1993). In fact, if one assumes 
that the anisotropy is constant throughout the system, the fourth 
order Jeans equation can be used to express the LOS velocity kur- 
tosis as an integral of the tracer density, anisotropy and total mass 
profiles (Lokas 2002). Moreover, Richardson & Fairbairn (2012) 
were able to generalize the expression for the LOS velocity kurto- 
sis for radially varying anisotropy in the framework of separable 
augmented density and 4th order anisotropy equal to the standard 
anisotropy (the latter appears to be an excellent approximation for 
ACDM haloes, see Fig. 10 of Wojtak et al. 2008). One can then per- 
form a joint fit of the observed LOS velocity dispersion and kurtosis 
profiles. This was found to (partially) lift the mass-anisotropy de- 
generacy when applied to dwarf spheroidal galaxies (Lokas 2002) 
and the Coma cluster, (Lokas & Mamon 2003): the joint constraint 
of LOS velocity dispersion and kurtosis profiles allows the estima- 
tion of both the mass profile (i.e., normalization and concentration) 
and the anisotropy of the cluster, contrary to the case when the LOS 
kurtosis profile is ignored. 

An interesting route is to perform non-parametric inversions 
of the data assuming either the mass profile to obtain the anisotropy 
profile (anisotropy inversion, pioneered by Binney & Mamon 1982) 
or the anisotropy profile to obtain the mass profile (mass inversion, 
independently developed by Mamon & Boue 2010 and Wolf et al. 
2010). These inversion methods are powerful in that they are non- 
parametric, but they suffer from their requiring the user to bin the 
data, smooth it, and extrapolate it beyond the range of data. 

Hence, one would like to go one step further and constrain the 
full information contained in the observed projected phase space 
(projected radii and LOS velocities, hereafter, PPS) of LOS veloc- 
ities as a function of projected radii. In other words, rather than 
using the 0th, 2nd and possibly 4th moments of the LOS velocity 
distribution, we wish to use the full set of even moments. 



The traditional way to analyze the distribution of particles in 
PPS is to assume a form for the six-dimensional distribution func- 
tion (DF) in terms of energy (E) and angular momentum (J) and 
fit the triple integral of equation (5) below, using this DF for /, 
to the distribution of particles in PPS. The worry is that we have 
no good a priori knowledge of the shape of the DF, f(E, J). One 
clever idea is to throw orbits in a gravitational potential, since each 
orbit is a Dirac delta function in energy and angular momentum. 
One then seeks a linear combination of these orbits, with positive 
coefficients, to match the data. This orbit model (Schwarzschild 
1979; Richstone & Tremaine 1984; Syer & Tremaine 1996) is very 
powerful (and can handle non-spherical gravitational potentials), 
but too slow to obtain meaningful errors on the parameters. A simi- 
lar, and in principle faster, technique is to assume that the DF is the 
linear combination (again with positive coefficients) of elementary 
DFs (Dejonghe 1989; Merritt & Saha 1993; Gerhard et al. 1998), 
but only one such study has been made (Kronawitter et al. 2000), 
and it is not clear whether the elementary DFs, although numerous, 
constitute a basis set. 

An important step forward has been performed by Wojtak 
et al. (2008), who analyzed the haloes in ACDM cosmological sim- 
ulations, to show that the DF can be approximated to be separable 
in energy and angular momentum, with a simple analytical approx- 
imation for the angular momentum term. In a sequel, Wojtak et al. 
(2009) have shown that it is feasible to fit the distribution of parti- 
cles in PPS with equation (5), using the approximation of the DF 
found by Wojtak et al. (2008). 

However, it is not yet clear whether self-gravitating quasi- 
spherical astrophysical systems have the DF of ACDM haloes: In 
particular, if the dynamical evolution of these systems is influenced 
by the dissipation of their gaseous component, the DF may not be 
separable in terms of energy and angular momentum. Dissipation 
is not expected to affect much the internal kinematics of large sys- 
tems such as galaxy clusters , but is expected to be increasingly 
important in smaller systems such as galaxy groups, and especially 
galaxies themselves. For this reason, it is useful to consider a mass- 
modeling method that is independent of the dependence of the DF 
on energy and angular momentum. 

In this work, we present an alternative method, in which we 
fit the distribution of particles in PPS making assumptions on the 
radial profiles of mass and anisotropy as well as the radial vari- 
ations of the distribution of space-velocities. We call this method 
Modelling of Anisotropy and Mass Profiles of Observed Spheri- 
cal Systems, or MAMPOSSt for short. 2 The MAMPOSSt method 
is described in Sect. 2.1, its Gaussian approximation is described 
in Sect. 2.2. Tests on haloes derived from a cosmological N body 
simulation are presented in Sect. 3. A discussion follows in Sect. 4. 



2 METHOD 

2.1 General method 

The observed tracer population of a spherical system has a DF 
f(r,v) =u(r)f v (v\r) , (2) 



1 However, the joint X-ray and lensing analysis of a cluster by Newman 
et al. (2009) reveals a shallower inner density profile than NFW, suggesting 
that dissipation is also important in clusters. 

2 MAMPOSSt should evoke the mass analog of a lamppost, and mam- 
posteria in Spanish means masonry, hence the building blocks of structures. 
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where v(r) is the tracer number density profile. 

MAMPOSSt fits the distribution of objects in PPS (projected 
radius R and LOS velocity v z ), assuming parametrized forms for 

(i) the gravitational potential (or equivalently a total mass den- 
sity profile, through the Poisson equation), 

(ii) the anisotropy profile (eq. [1]) 

(iii) the distribution of 3D velocities, f v (v\r). 

Consider a point P at distance r from the centre, O, of the 
spherical system, with projected radius R ^ r and consider the 
spherical coordinates where the unit vectors e r and eg are in the 
plane containing OP and the LOS, while is perpendicular to this 
plane. Consider also the cylindrical coordinate system (v z , v± , v^), 
where e, is the axis along the LOS and e± is the axis perpendic- 
ular to the LOS, but in the plane containing O and P and the LOS. 
The Jacobian of the transformation from the spherical coordinate 
system to the new one is unity, hence one can write 



fv(v z ,v±,v lj> \{r,R}) = 



d d N 



dv z dv± di>0 



d d N 



dv r dvg di)0 



= fv(v r ,v e ,v^\{r,R}) . 

The distribution of LOS velocities at P is then obtained by 
integrating velocities over the two perpendicular axes (dropping 
{r, R} from /„ for clarity): 



h(v z \R,r) 



' dN 





r+oo 




^ r,R 


/ dv ± 


I 




1 —CO 





dv_L / fv(v z , v±,v <t> ) dv,f, .(3) 



Note that dynamical systems have maximum velocities set by the 
escape velocity, — 2 <£>(r) (where $(r) is the gravitational po- 
tential), on one hand, and by the maximum allowed (observable) 
absolute LOS velocity on the other hand. In what follows, we will 
neglect both limits, unless explicitly mentioned otherwise. 

The surface density of observed objects {the tracer) in PPS is 
then obtained by integrating along the LOS 



g(R,v z ) = E(R) (h(v z \R,r)) hOS 
r v(r) 



r Vr 



R 2 



h(v z \R,r) dr 



(4) 



r dr 



where 

E(fl) = 



R \fr T ~^R? 



v(r)dz = 2 



dv± f(r, v z ,v±, v^dv^ .(5) 



(6) 



r v(r) dr 



R 2 



is the tracer surface density at projected radius R. Equation (5) is 
equivalent to equation (2) of Dejonghe & Merritt (1992). 

If the tracer number density profile v(r), appearing in equa- 
tion (4), is not known and if the incompleteness of the data is in- 
dependent of the projected radius, then one can estimate u(r) by 
Abel inversion of E(ii) of equation (6): 



v{r) = -- 



dE dR 



7T J r dR 



(7) 



But this is not necessary, as we shall see below. 

In MAMPOSSt, rather than replace the velocities by energy 
and angular momentum and numerically solve the triple integral 
of equation (5) (as first proposed by Dejonghe & Merritt, see also 



Wojtak et al. 2009), we analytically derive h(v z \R, r) from equa- 
tion (3) for known 3D velocity distributions. With the analytical 
form of h(v z \R, r), equation (4) provides the surface density dis- 
tribution of tracers in PPS through a single integral. Note, however, 
that another single integral is required because the expression for 
h(v z \R, r) will involve ay (r) (see eqs. [25] and [26], below, for the 
Gaussian case), which is obtained by solving the spherical Jeans 
equation 



dr 



+ 2/3- 



-u{r) 



GM(r) 



(8) 



where /3 is the anisotropy of (eq. [1]) for our given choices of total 
mass and anisotropy profiles. We thus need to insert the solution 
(van der Marel 1994; Mamon & Lokas 2005) 



(r) 



1 

v(r) 



exp 



t s <fM{s) 
v(s) 5 — ds , 



(9) 



in the expression for h(v z \R, r) (eq. [3]) to derive g(R,v z ), via 
equation (4), where /3(f) is given, while u(r) is obtained with equa- 
tion (7). In equation (9), M (s) = (s 2 /G) d$/ds is the radial pro- 
file of the total mass (this is the only instance where the gravita- 
tional potential enters MAMPOSSt). For a given choice of parame- 
ters, the single integral of equation (4) must be evaluated for every 
data point (R,v z ), whereas the other integral (eq. [9]) for oy(r) 
need only be evaluated once, on an adequate grid of r. 

Note that for projected radii extending from R m in to 7? ma x 
and absolute LOS velocities extending from to a maximum veloc- 
ity, which for projected radius 7? is theoretically equal to v csc (R) = 
y/—2<f>(R), and in practice is possibly specified by a cut of obvi- 
ous velocity interlopers, v cut (R), one can write 



2nRdR 



g(R,v z ) dv z 

V cu t-R) 



2tt 



A7V D 



RE(R) dR 
(10) 



where we used equation (4) for g(R, v z ), assumed that h(v z \ R, r) 
is normalised, reversed the order of the integrals in r and v z , and 
where N P (R) is the predicted number of objects within projected 
radius R, while AN P = 7Yp(7? ma x) - A r P (i?min). Equation (10) 
then implies that the probability density of observing an object at 
position (R, v z ) of PPS is 

2nRg(R,v z ) 



q(R,v z 



AN D 



4nR 
AiVp 

R 2 
A7V p rg 



r v(r) 



oo 

cosh u v 



R- 
R 



h(v z \R,r)dr (11) 

cosh iij h(v z \R, i?cosh u) du, 
(12) 



where equation (11) arises from equation (4), while equation (12) 
is obtained by writing r = i?coshti. Here, N P (R) is the number 
of tracers in a cylinder of projected radius R, the terms v and 7V P 
are given by 



u(r) = N ^ v( r 
Aivrl \r„ 

N p (R) = N(r„) N p ( — 



(13) 
(14) 



where N(r) is the cumulative tracer number density profile, while 
r v is the characteristic radius of the tracer. One easily verifies that 
J J q(R, v z ) dRdv z — 1. The values of R m in and i? max appearing 
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in ATV P (eq. [10]) can be hard limits, or alternatively the respective 
minimum and maximum projected radii of the observed tracers if 
no hard limits are specified. 

We fit the parameters (mass scale or concentration and possi- 
bly normalization, anisotropy level or radius, as well as the tracer 
scale r v - if not previously known) that enter the determination of 
g(R,v z ) to the observed surface density, using maximum likeli- 
hood estimation (MLE), i.e. by minimizing 



ln£ = -^>2lnq(Ri,v z ,i\0) 



(15) 



for the TV-parameter vector 9, where n is the number of data points, 
with q given by equation (11). 

Writing — {r„,Tj}, where tj is the vector of the TV — 1 
parameters other than r v , one has 



q(R,v z \r v , tj) = p (R\r v ) x p{v z \R,r v ,rj) , 



where 



po(R\r„) 



and 



2tt RE(R\r„) 



N p (R n 



p(v z \R,ri) = (h(v z \R,r)) 1 



N p (Rmin\r v ) 



9(R,v z \0) 



(16) 



(17) 



(18) 



Combining the last equality of equation (18) with equation (5), in- 
tegrating over LOS velocities, reversing the order of the two outer 
integrals of the resulting quadruple integral, and using equation (6) 
yields J p(v z \R) dv z — 1. So, if the scale of the tracer distribu- 
tion is already known, then, according to equations (15) and (16), 
maximizing the likelihood amounts to minimizing 



ln£' = — 2^ lnp(^z|-R, rv, tj) 



(19) 



Now, if r„ is not known, then one may be tempted to solve for 
it by minimizing — ln£o = — m Po(-R|Tv), an d then proceed 
with equation (19) to minimize for the TV — 1 remaining parame- 
ters, tj. However, since — In C = — In £' — In Co (from eqs. [15] 
and [16]), the most likely solution for that minimizes — In C will 
not in general be that which minimizes at the same time — ln£' 
and — lnXn- Moreover, if one seeks to obtain the distributions of 
parameters tj and r„ consistent with the MLE solution (for exam- 
ple with Markov Chain Monte-Carlo techniques), the joint analysis 
of equations (11) and (15) is required. On the other hand, if r v is 
known from other data, while the current dataset is known to have 
a completeness, C(R), that is a function of projected radius, then 
one could indeed minimize ln£' of equation (19). The proper so- 
lution is then to minimize — In £ weighting the data points by the 
inverse completeness, i.e. minimizing 



InC" 



E 



In q(Rj, v z ,i\0) 
C(Ri) 



(20) 



For computational efficiency, we perform the following tasks: 

(i) For each run of parameters, we first compute log ay (rj ) from 
equation (9) on a logarithmic grid of Tj, and compute cubic-spline 
coefficients at these radii. Then, when we compute the LOS in- 
tegral of equation (4) for each (Rj,v z ,i), we evaluate a r (r) with 
cubic spline interpolation (in log-log space, using the cubic spline 
coefficients determined at the start). 

(ii) For simple anisotropy models, the exponential term in equa- 
tion (9) is given by equations (A2) and (A3). 



(iii) We terminate the LOS integration in equation (11) at 
roughly 15 virial radii, 3 r v , instead of infinity, as the Hubble flow 
pushes the velocities of the material beyond this distance to val- 
ues over 3 a v above the mean of the system (see Mamon, Biviano, 
& Murante 2010, hereafter MBM10). The LOS integration varies 
only very slightly with the number of virial radii, so as long as the 
virial radius is correct to a factor of two, this choice of integration 
limit is not an issue. 

We now need to choose a model for the shape of the 3D veloc- 
ity distribution. While MAMPOSSt, can, in principle, be run with 
any model, the simplest one is the (possibly anisotropic) Gaussian 
distribution, which we describe in Sect. 2.2 below. 



2.2 Gaussian 3D velocity distributions 

The simplest assumption for the 3D velocity distribution is that it 
is Gaussian: 



fv(Vr,V ,V^) = 



(27r) 3 / 2 a r <jj 



■ exp 



2cr 2 



4 + tj 



(21) 



where the velocity dispersions cr, are functions of r. This Gaussian 
distribution assumes no streaming motions: e.g. no rotation, and 
no mean radial streaming, which is adequate for i? max < r v in 
high-mass haloes (i.e. groups and clusters) and i? max < 4r v in 
galaxy-mass haloes (Cuesta et al. 2008). Inserting equation (21) 
into equation (3) and integrating over v$ leads to 



h(v z \R,r) = 



, 2-Ky/T~P a; 



exp ■ 



V 2 r +Vl 

2(1-/3) a? 



(22) 



Calling 9 the angle between the line-of-sight (direction z) and the 
radial vector r, one has 



v r = v z cos + v± sin , 
v$ = —v z s'm8 + v± cos ( 



(23) 
(24) 



with which the integral over v± in equation (22) yields a Gaussian 
distribution of LOS velocities at point P: 

h(v z \R, r) = — : — exp 



of squared dispersion 



2al(R,r) 



(25) 



cr 2 (R,r) = 



l-«r)(|) 



a 2 (r) 



(26) 



The integral of h(v z \R,r) along the LOS is obtained from equa- 
tions (4) and (25): 



g(R,v z 



l-/3R 2 /r 2 



7T J R V^ 2 - R? 

x exp 



2 (l-/3R 2 /r 2 ) al 



dr . (27) 



According to equations (18) and (27), the probability of measuring 
a velocity v z at given projected radius R is 



p(v z \R) = 



E(fl) 



3 The virial radii are loosely denned here as the radius where the mean 
density of the halo is 200 times the critical density of the Universe. 
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f™(y/<r,) exp [-v 2 J (2ol)\ dz 



>2k 



r 

Jo 



v dz 



(28) 



We remind the reader that /3 is a chosen function of r, v is 
a function of r given by equation (7), while o r is a function of 
r given by equation (9). For isotropic systems (/? = 0), equa- 
tion (27) leads to a Gaussian distribution of LOS velocities. How- 
ever, for anisotropic velocity tensors, the distribution of LOS ve- 
locities will generally not be Gaussian (as Merritt 1987 found when 
starting from distribution functions instead of Gaussian 3D veloc- 
ities). Hence, the Gaussian nature of h(v z \R, r) is not equivalent 
to the popular assumption that g(R, v z ) is Gaussian on v z : even if 
h(v z \R, r) is a Gaussian at point P, its integral along the LOS is 
not Gaussian, unless /3 — and ay is constant. 

If one of the parameters to determine with MAMPOSSt is the 
normalization of the mass profile, one should not be tempted in 
expressing the radii in terms of the virial radius r v , the velocities in 
terms of the virial velocity v v , the tracer densities in terms of what 
we wish (as they appear in both the numerator and denominator of 
eq. [28]). Doing so, equation (28) becomes 



p(vz\R) 



v v p(v z \R) 
1 f °° ex P dz 



/2t? 



JT udz 



(29) 



where the quantities with tildes are in virial units. Equation (29) 
indicates that when one varies the r v (and the virial velocity in 
proportion as v v = -y/2/A Ho r v ), the highest probabilities are 
reached for the highest normalizations: v z becomes very small, 
while o~ z is unaffected to first order. This unphysical result is the 
consequence of using a parameter (the virial radius) as part of the 
data variable. On the other hand, using equation (28), one sees that 
the highest probabilities p(v z | R) are reached at intermediate values 
of the normalization. 

Taking the second moment of the velocity distribution of equa- 
tion (28) leads to the equation of anisotropic projection yielding the 
LOS velocity dispersion, a z (R): 



T,(R)a 2 (R) = 



\g{R,v z ) di 



vr dr 



n Jr <r z (R,r) Vr 2 - R 2 



v 2 exp 



2a 2 (R,r) 



dv z 



2 / 

' R 



0(r) 



R' 2 



r dr 



R? 



(30) 



Equation (30) recovers the equation of anisotropic kinematic pro- 
jection, first derived by Binney & Mamon (1982). 

If interlopers are removed with a velocity cut v cu t(R), then 
the expression for h(v z \R, r) becomes 



h(v z \R,r) = 



exp {-v 2 / [2<T 2 (R,r)}} 



2na z (R,r)erf{v cut (R)/ [a z (R, r)y/2~\ } 



•(31) 



In summary, MAMPOSST with Gaussian 3D velocities com- 
putes likelihoods from equations (15), (11) or (12), (25), (26), and 
(9), in that order. 



3 TESTS 

3.1 Simulated haloes 

To test MAMPOSSt, we use cluster-mass haloes extracted by Bor- 
gani et al. (2004) from their large cosmological hydrodynamical 
simulation performed using the parallel Tree+SPH GADGET-2 code 
of Springel et al. (2005). The simulation assumes a cosmological 
model with fi = 0.3, Ov = 0.7, f2 b = 0.039, h = 0.7, and 
erg = 0.8. The box size is L = 192 h^ 1 Mpc. The simulation used 
480 3 DM particles and (initially) as many gas particles, for a DM 
particle mass of 4.62 x 10 9 /i _1 M0. The softening length was set 
to 22.5 h^ 1 comoving kpc until z = 2 and fixed afterwards (i.e., 
7.5 hT 1 kpc). The simulation code includes explicit energy and 
entropy conservation, radiative cooling, a uniform time-dependent 
UV background (Haardt & Madau 1996), the self-regulated hy- 
brid multi-phase model for star formation (Springel & Hernquist 
2003), and a phenomenological model for galactic winds powered 
by Type-II supernovae. 

DM haloes were identified by Borgani et al. (2004) at redshift 
z = with a standard Friends-of-friends (FoF) analysis applied 
to the DM particle set, with linking length 0.15 times the mean 
inter-particle distance. After the FoF identification, the centre of the 
halo was set to the position of its most bound particle. A spherical 
overdensity criterion was then applied to determine, for each halo, 
our proxy for the virial radius, r2oo, where the mean density is 200 
times the critical density of the Universe. 

To save computing time, we worked on a random subsample 
of roughly two million particles among the 480' ! . We have extracted 
1 1 cluster-mass haloes from these simulations, among which, ten 
are about logarithmically spaced in virial radius, r2oo, while the 
11 th halo is the most massive in the entire simulation. Their prop- 
erties are listed in Table 1. We made no effort to omit irregular 
haloes, but among the list of 12 irregular haloes out of 105 ex- 
tracted by MBM10 from the same simulation, 2 are in our sample 
(haloes 17283 and 434). We list the characteristic radii r s , ru, tb 
of three models fitted by MLE to the mass density profiles of the 
particle data (from 0.03 to 1 r2oo), namely: 



(i) the NFW density profile 

p{r) oc r^ 1 (r + r s )~ 2 , 



(32) 



where r s = r_2 is the radius of slope —2 in the mass density 
profile, related to the concentration c = r2oo A-2 ; 
(ii) the Hernquist density profile (Hernquist 1990) 



p(r) oc r 1 (r + r H ) ' ! , 



where ru — 2 r_2, 

(iii) the Burkert density profile (Burkert 1995) 

p(r) oc (r + r B ) _1 (r 2 + rl) 1 , 
where tb — 0.657 r_2. 



(33) 



(34) 



Denoting the scales r s , ru and re by the generic r p , the mass pro- 
files of these models (required for eq. [9]) can be written 



M(r) = M{r p ) 



where 



M (r/r p 
M(l) 



(35) 
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Table 1. Properties of 1 1 cosmological haloes 



rank 


ID 






? 'T T 


i'T> 


A 


To 


Aoo 


1 


18667 


0.789 


0. 179 


0.401 


0.117 


1.14 


U.Z/6 


1 .33 


2 


21926 


0.842 


0.123 


0.342 


0.085 


1.34 


0.050 


1.73 


3 


30579 


0.890 


0.189 


0.443 


0.120 


1.33 


0.053 


1.69 


4 


25174 


0.956 


0.144 


0.377 


0.099 


1.23 


0.162 


1.36 


5 


3106 


1.010 


0.297 


0.661 


0.166 


1.05 


2.384 


1.09 


6 


8366 


1.076 


0.434 


0.819 


0.249 


1.11 


0.689 


1.29 


7 


13647 


1.151 


0.227 


0.536 


0.151 


1.19 


0.265 


1.41 


8 


1131 


1.174 


0.197 


0.499 


0.133 


1.18 


0.352 


1.34 


9 


17283 


1.298 


0.505 


1.009 


0.277 


1.04 


0.727 


1.05 


10 


434 


1.374 


0.317 


0.699 


0.210 


1.30 


0.165 


1.70 


1 1 


5726 


1.660 


0.407 


0.921 


0.249 


1.42 


0.050 


2.20 



Stack 1.09±0.08 0.26±0.04 0.60±0.08 0.17±0.02 1.21±0.04 0.26±0.08 1.45±0.10 



Notes: Properties obtained from fits to the particle data of 1 1 haloes. Cols. 1 and 2: cluster identification; col. 3: virial radius T200\ col. 4: scale radius (= r_2) 
of the NFW mass density profile (eq. [32]); col. 5: scale radius (= 2 r_2) of the Hernquist mass density profile (eq. [33]); col. 6: scale radius (~ 0.657 r_2) 
of the Burkert mass density profile (eq. [34]); col. 7: mean anisotropy (.A = cr r /crg) within r200\ col. 8: anisotropy radius with the ML anisotropy model; 
col. 9: asymptotic anisotropy (Aoo = crr/crg) at infinite radius with the T anisotropy model. Radii are in units of h~ 1 Mpc. The measured anisotropies do 
not incorporate streaming motions. 



'ln(a; + l)- 



x + 1 



(NFW) 



M(x) = 



\x + lj 



(Hernquist) ( 36 ) 



M [(a; + l) 2 (a; 2 + l) 



2 tan x . (Burkert) 



The NFW model has long been known to fit well the den- 
sity profiles of ACDM haloes (Navarro et al. 1996), and while 
Navarro et al. (2004) found that Einasto models fit them even bet- 
ter, MBM10 found that the NFW model describes the outer LOS 
velocity dispersion profile of the DM component of their stacked 
cluster-mass halo in Borgani et al.'s hydrodynamical cosmologi- 
cal simulation even (slightly) better than the Einasto model. The 
Hernquist model differs from the NFW one because it has a steeper 
logarithmic slope at large radii, 7 = d In p/d In r — —4 rather than 
—3. The Burkert model, on the other hand, has the same asymptotic 
7 = —3 as the NFW model, but a core at the centre, 7 = 0, rather 
than a cusp (7 = — 1 in both the NFW and Hernquist models). 

In Table 1 , we also list the values of the parameters character- 
izing different velocity-anisotropy models, namely: 

(i) the constant anisotropy model a r /ag = (1 — /J) -1 / 2 = A 
('Cst' model hereafter), where we assume spherical symmetry and 
therefore erg = a$; 

(ii) the model ('ML' model hereafter) of Mamon & Lokas 
(2005); 



1 r 
/3ml (r) = - — ■ 

2 r + ra 



(37) 



characterized by the anisotropy radius rp\ 

(iii) a generalization of the ML model, which is also a simplified 
version of the model of Tiret et al. (2007), isotropic at r — and 
with anisotropy radius identical to r_2 (hereafter called 'T' model): 



Pr(r) = Po 



r + r-2 



(38) 



characterized by the anisotropy value at large radii, fioo- In our T 
model, the anisotropy radius is set to the scale radius of the mass 
density profile. Note also that in the following we provide the val- 
ues of Aoo = (cr r /rjg) 00 = (1 — /3oo)~ 1/2 , rather than /?oo. 




Figure 1. Velocity anisotropy profiles of the 1 1 haloes (broken coloured 
lines). The smooth black curve is the ML anisotropy model with rp = r-2 
(or, equivalently, the T anisotropy model with 0^ = 0.5). 



Note that the ML and the T models used here are identical for 
Poo =0.5 and rp = r_2- With these values, the ML and T mod- 
els provide a good fit to the average anisotropy profile of a set of 
cluster-mass cosmological haloes (MBM10). 

In Fig. 1, we show the individual halo velocity anisotropy pro- 
files and the ML anisotropy model with r$ = r_2 (or, equivalently, 
the T anisotropy model with j3oo — 0.5). There is a huge scatter in 
the P(r) of the individual haloes, as already observed by, e.g., Wo- 
jtak et al. (2008), especially at r > 0.3r2oo, while 8 out of 11 
haloes have ,3(0) = ± 0.15. 



© 2012 RAS, MNRAS 000, 000-000 



Modelling Anisotropy and Mass Profiles 7 



3.2 Observing cones and interloper removal 

To test MAMPOSSt, we select 500 DM particles around each halo, 
out to a maximum projected distance i? m ax from the halo centre, 
for which we consider three values: rsoo — 0.66 V2oo, ^200, and 
rioo — 1-35 r2oo- We analyze three orthogonal projections for each 
halo - these are in fact cones with an observer at D — 90 h^ 1 Mpc 
away, but the opening angle being very small has no noticeable ef- 
fect on our results. The particles in these cones are used by MAM- 
POSSt as tracers of the halo gravitational potential. 

However, these 500-particle samples include interlopers, i.e. 
DM particles that are located in projection at R $C i? m ax, but are 
effectively outside i? max in real (3D) space, i.e. with r > i? ma x- It 
is impossible to remove all these interlopers in the observed redshift 
space, where only 3 of the 6 phase-space coordinates of the tracers 
are known (e.g., MBM10). Moreover, since the LOS velocity dis- 
tribution of interlopers in mock cones around ACDM haloes is the 
sum of a Gaussian component and a uniform one (see MBM10 for 
a quantified view), and since the Gaussian one resembles that of the 
particles in the virial sphere, it is important to remove the flat LOS 
velocity component, at least at high absolute LOS velocity, where 
it dominates. It is possible to remove these high \v z \ objects with 
suitable interloper removal algorithms. 

To see how MAMPOSSt depends on the choice of the in- 
terloper removal algorithm, we here consider three different algo- 
rithms. 

The first one is a new, iterative algorithm, that we name 
"Clean", which is fully described in Appendix B. Clean first looks 
for gaps in the LOS velocities, then estimates the virial radius V2oo 
from the aperture velocity dispersion, assuming an NFW model 
with ML anisotropy with an anisotropy radius rp = r_2 and a 
concentration depending on the estimate of r2oo via the relation of 
Maccio, Dutton, & van den Bosch (2008), then only considers the 
galaxies within 2.7cr z (R) from the median LOS velocity, and fi- 
nally iterates. Our assumed anisotropy profile fits reasonably well 
the anisotropy profiles of DM haloes (MBM10), as is clear for our 
1 1 haloes (see Fig. 1). The factor 2.7 was found by MBM10 to best 
preserve the local LOS velocity dispersion for the assumed density 
and anisotropy models. 

We also consider two other interloper removal algorithms, 
namely: 

(i) the method (hereafter, dHK) of den Hartog & Katgert (1996), 
a widely used procedure that works reasonably well on cluster- 
mass haloes from cosmological simulations (Biviano et al. 2006; 
Wojtak et al. 2007), despite its crude underlying physics; 

(ii) the method (hereafter, KBM) of Katgert, Biviano, & Mazure 
(2004, see their Appendix A), in which a galaxy is flagged as an 
interloper under the condition v z /a z > 1.85 (-R/ r 20o)~ 0,3 , with 
r2oo derived from a z using eq. (8) of Carlberg, Yee, & Ellingson 
(1997). This method was invented as a poor-man's proxy for the 
dHK method when the observational sampling of the halo projected 
phase-space is poor. 



3.3 The general 4-parameter case 

There is no a priori limitation on the number of free parameters 
that can be used in MAMPOSSt to characterise the mass and ve- 
locity anisotropy profiles. With samples of ^ 500 tracers (assumed 
massless throughout these tests) it is appropriate to consider ~ 4 
free parameters, two for the mass distribution, one for the velocity 
anisotropy distribution, and one for the spatial distribution of the 



tracers. All these models are characterised by the two free parame- 
ters, the 'virial' radius r2oo and a characteristic scale-radius (r s , th, 
and re for the NFW, Hernquist, and Burkert models, respectively). 
Herafter, we generically use r p to refer to this characteristic scale- 
radius of the mass density profile. 

We use the NFW model, in projection (Bartelmann 1996; 
Lokas & Mamon 2001), to fit the projected number density pro- 
file of the tracer. Note that the normalization of this profile does 
not enter the MAMPOSSt equations, so the only free parameter is 
r_2. Herafter we call this parameter r v , to avoid confusion with 
the characteristic radius of the NFW mass density profile. We only 
consider one model for the number density profile of the tracer, be- 
cause this is a direct observable, unlike the mass density profile. 
While one should not be too restrictive in the model choice for the 
mass density profile, the observer is generally able to choose the 
best-suited model for the tracer number density profile by direct 
examination of the data before running MAMPOSSt. We choose 
the NFW model because it provides a reasonable description of 
the number density profiles of the DM particles in our simulated 
haloes. 

For the velocity anisotropy profile, we consider the three mod- 
els described above, Cst, ML, and T, each characterised by a sin- 
gle anisotropy parameter, A,rp, and A x , respectively. In equa- 
tion (38), we use r_2 = r p . 

To search for the best-fit solution, we run the MAMPOSSt 
algorithm in combination with the NEWUOA 4 minimization rou- 
tine of Powell (2006). For estimating error bars on the best fit pa- 
rameters, as well as confidence contours on pairs of parameters, 
we fit our model parameters using the Markov Chain Monte Carlo 
(MCMC) technique (e.g., Lewis & Bridle 2002). In MCMC, the k- 
dimensional parameter space is populated with proposals, for each 
of which the likelihood is computed. The new proposal is accepted 
if the ratio of new to previous likelihood is either greater than unity 
or else greater than a uniform [0, 1] random number. The proposal 
is found by assuming a fc-dimensional Gaussian probability distri- 
bution around the previous proposal. We adopt the publicly avail- 
able CosmoMC code by A. Lewis. 5 We run 6 chains in parallel 
using Message Parsing Interface (MPI), and the covariance matrix 
is used to update the parameters of the Gaussian proposal density 
to ensure faster convergence. 

Fig. 2 illustrates the MAMPOSSt analysis via MCMC for 
the general case with four free parameters, using the NFW model 
for the mass density profile, and the constant (free parameter) 
anisotropy model. In particular, it shows that the different parame- 
ters are not correlated, except for a positive correlation between r p 
and A. 

Our FORTRAN code takes roughly 1 second to find the 
MAMPOSSt 4-parameter solution for a 500 particle sample run 
in scalar on a decent desktop or laptop computer, and 4 minutes 
to produce confidence limits for this solution with the CosmoMC 
(Lewis & Bridle 2002) MCMC code, with 6 chains of 40 000 ele- 
ments run in parallel (MPI) on a PC equipped with a 4-core 8-thread 
Intel Core-I7 2600 processor. 

The results for the different interloper rejection methods, mass 
density and velocity anisotropy models, and for the different max- 
imum projected radii used in the selection of the 500 particles are 
listed in Table 2 and displayed in Fig. 3. We list and show the 
biweight measures (see, e.g., Beers, Flynn, & Gebhardt 1990) of 



4 NEWUOA is available at http://plato.asu.edu/ftp/other_software/newuoa.zip 

5 CosmoMC is available at http : / /co sinologist . inf o/cosmomc/. 
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Figure 2. Illustration of MAMPOSSt analysis for the general case, with ft independent of radius, for a 500 particle sample from axis x of halo 25174 (grey 
broken line in Fig. 1), using MCMC (with 6 chains of 40 000 elements). The contours are 1, 2, and 3cr. The red arrows and stars indicate the maximum 
likelihood solution, while the green arrows and crosses show the true solution (Table 1). The priors for the MCMC were uniform within the boxes of each 
panel and zero beyond the boxes. 



mean and dispersion of log(o/t) where o is the recovered value of 
the parameter and t its true value, because, according to our tests, 
they perform better than standard statistical estimators of location 
and scale when the parent distributions are not pure Gaussians. We 
call 'bias' and 'inefficiency' the mean and dispersion of \og(o/t). 
If the dispersion in true values of a given parameter is small, one 
can spuriously obtain low values of the log(o/t) dispersion when 
the MAMPOSSt and true values show no correlation. We therefore 
also list the Spearman rank correlation coefficient between o and 
t, marking in boldface those correlations that are significant at the 
99% confidence level. We list the results for all interloper rejection 
methods and all anisotropy models only for the NFW mass density 
model and for the R ^ V2oo radial selection. For simplicity, we 
only show a limited set of results for the other mass density models 
and for the other radial selections. 

Remarkably, as seen in Table 2, the results for the four param- 
eters are almost independent of the interloper removal algorithm, 
the Clean and KBM algorithm performing slightly better than dHK. 
The results for T2oo , r v , and r p also depend very little on the cho- 



sen anisotropy model. On average, the values of the raoo parameter 
are recovered with almost no bias (from —1 to +4%) and with only 
~ 10% inefficiency. The r v parameter estimates are always slightly 
positively biased (4—7%), and are recovered with ~ 25% efficiency. 
Also, the r p parameter estimates generally display a slight positive 
bias, except for the KBM interloper removal method, and overall 
the bias ranges from —2 to +15%, while the efficiency ranges from 
~ 50 to ~ 90%. 

As far as the anisotropy parameter is concerned, the ML model 
behaves very differently from the Cst and T models, in that it is 
virtually impossible to constrain the anisotropy radius of the for- 
mer, rp, while it is possible to obtain quite good constraints on 
the anisotropy parameters of the other two models, A and Aoa- 
More precisely, the rp estimates are always negatively biased (by 
~ 60%) and are affected by a huge dispersion (almost one order of 
magnitude). On the other hand, A and Aac are almost unbiased (the 
bias ranges from —10 to +5%) and they are affected by dispersions 
of, typically, ~ 20%, if we consider the Clean and KBM interloper 
removal algorithms. So, apparently, it is much easier to constrain 
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Table 2. MAMPOSSt results for different interloper removal algorithms, density models, apertures, and number of particles 



N 


^max 


Membership 


p(r) 


0(r) 




T200 






r v 






r p 




anisotropy 














bias 


ineff. 


COIT. 


bias 


ineff. 


corr. 


bias 


ineff. 


corr. 


bias 


ineff. 


corr. 


500 


f200 


Clean 


NFW 


Cm 


0.004 


0.040 


0.909 


0.027 


0.102 


0.835 


0.032 


0.217 


0.578 


0.007 


0.073 ■ 


-0.255 


500 


^200 


Clean 


NFW 


ML 


-0.003 


0.040 


0.904 


0.024 


0.104 


0.832 


0.057 


0.229 


0.601 


-0.221 


0.887 


-0.172 


500 


^200 


Clean 


NFW 


T 


-0.006 


0.040 


0.903 


0.026 


0.103 


0.838 


0.039 


0.169 


0.709 


0.007 


0.085 


0.621 


500 


^200 


dHK 


NFW 


Cst 


0.018 


0.042 


0.885 


0.027 


0.099 


0.838 


0.051 


0.319 


0.406 


0.004 


0.147 


-0.215 


500 


r200 


dHK 


NFW 


ML 


0.012 


0.041 


0.909 


0.028 


0.100 


0.840 


0.059 


0.286 


0.611 


-0.161 


0.904 


-0.1 18 


500 


''200 


dHK 


NFW 


T 


0.012 


0.044 


0.902 


0.027 


0.100 


0.844 


0.022 


0.199 


0.636 


-0.045 


0.264 


0.464 


500 


^200 


KBM 


NFW 


Cst 


0.005 


0.038 


0.906 


0.018 


0.100 


0.850 


-0.006 


0.218 


0.535 


0.020 


0.078 


-0.198 


500 


^"200 


KBM 


NFW 


ML 


-0.003 


0.039 


0.908 


0.020 


0.100 


0.851 


-0.005 


0.232 


0.557 


-0.191 


0.795 


0.101 


500 


^200 


KBM 


NFW 


T 


-0.006 


0.038 


0.911 


0.020 


0.099 


0.856 


-0.010 


0.184 


0.689 


0.018 


0.094 


0.595 


500 


^200 


Clean 


Her 


T 


0.002 


0.039 


0.909 


0.026 


0.102 


0.835 


0.039 


0.132 


0.755 


0.014 


0.086 


0.546 


500 


^200 


Clean 


Bur 


T 


0.000 


0.039 


0.910 


0.047 


0.196 


0.377 


0.048 


0.145 


0.704 


-0.019 


0.071 


0.603 


500 


T500 


Clean 


NFW 


T 


-0.004 


0.048 


0.877 


0.089 


0.115 


0.902 


0.016 


0.143 


0.744 


0.009 


0.108 


0.232 


500 


^100 


Clean 


NFW 


T 


-0.014 


0.035 


0.905 


0.039 


0.179 


0.420 


0.093 


0.210 


0.538 


-0.001 


0.090 


0.436 


100 


T200 


Clean 


NFW 


Cst 


-0.001 


0.058 


0.844 


0.033 


0.201 


0.537 


-0.133 


0.341 


0.342 


0.003 


0.119 


-0.053 


LOO 


^200 


Clean 


NFW 


ML 


-0.011 


0.061 


0.834 


0.034 


0.199 


0.539 


-0.087 


0.336 


0.484 


-0.137 


0.925 


-0.245 


LOO 


^200 


Clean 


NFW 


T 


-0.008 


0.058 


0.850 


0.032 


0.200 


0.532 


-0.108 


0.277 


0.436 


0.014 


0.143 


0.249 



Notes: These results are for 1 1 haloes each observed along 3 axes, general 4 free-parameter case. Col. 1 : Number of initially selected particles (before interloper 
removal); col. 2: Maximum projected radius for the selection, where rsoo — 0.65 r2uo and ^100 — 1.35r20o; col. 3: Interloper-removal method (dHK: den 
Hartog & Katgert 1996; KBM: Katgert et al. 2004; Clean: App. B); col. 4: mass density model (NFW: Navarro et al. 1996; Her: Hernquist 1990; Bur: Burkert 
1995); col. 5: anisotropy model (Cst: (3 = cst; ML: eq. [37], Mamon & Lokas 2005; T: eq. [38], adapted from Tiret et al. 2007); cols. 6-8: virial radius; 
cols. 9—11: tracer scale radius; cols. 12—14: dark matter scale radius; cols. 15-17: velocity anisotropy (i.e., A for the Cst model, rp for the ML model, and 
Aoa for the T model). The columns 'bias' and 'ineff.' respectively provide the mean and standard deviation (both computed with the biweight technique) of 
log(o/i), while columns 'corr.' list the Spearman rank correlation coefficients between the true values and MAMPOSSt-recovered ones (values in boldface 
indicate significant correlations between o and i values at the > 0.99 confidence level). 



the 'normalisation' of a given anisotropy profile, than to constrain 
the characteristic radius at which the anisotropy changes, particu- 
larly so if this change is mild, as in the ML model. Note, however, 
that the difficulty of MAMPOSSt in constraining the anisotropy 
parameter of the ML model does not mean that the ML model is 
a poor representation of reality, and in fact Fig. 1 suggests the op- 
posite. Moreover, constraints obtained on the r200,?V, and r p pa- 
rameters are equally good with the Cst and T anisotropy models, as 
with the ML one. 

As seen in Table 2, correlations between recovered and ob- 
served values of the parameters r2oo, r v , and r p are almost always 
signficant. This is also true for the Aoo anisotropy parameter, but 
not for A and 773. In Fig. 4, we show the correlations existing be- 
tween the true and recovered values of the different parameters, us- 
ing the T anisotropy model, for the 1 1 haloes along the 3 different 
orthogonal projections. Projections effects render the determination 
of the mass and anisotropy profile of a single 500-particle halo very 
uncertain. However, Fig. 4 shows that ~ 500 tracers are sufficient 
to rank haloes for each of the different parameters considered here, 
i.e. by mass 0*200), scale radius of the tracer distribution (r v ) and of 
the total mass distribution (rp), and outer velocity anisotropy Aoa- 

The importance of projection effects is also very clear from 
Fig. 5, where we display the ratio of the recovered to true values 
of the parameters for each halo along the three different projection 
axes. This figure also shows that there is no trend of under- or over- 
predicting the parameter values with halo mass. 

All the above considerations apply for the NFW mass pro- 
file. Our tests with the Hernquist and Burkert mass profiles give 
similar results, as can be seen in Table 2, where for simplicity, we 
only list the results for the Clean interloper-removal algorithm and 



for the Tiret anisotropy model. The results are very similar for the 
NFW and Hernquist models. Results are similar also for the Burk- 
ert model, except for the scale r„ of the number density profile, 
for which the bias and inefficiency are both higher than those ob- 
tained using the NFW and Hernquist mass models. Since the model 
we use for the number density profile has not changed (a projected 
NFW), this result suggests that it is difficult to accomodate a tracer 
with a central cuspy spatial distribution in a potential with a central 
core. 

All the results described so far were obtained for a selection 
of 500 particles within i? max = ?*200 • Changing the value of i? max 
is not without effects on the results. The inefficiency on r2oo de- 
creases when i? m ax is gradually increased from rsoo to noo- The 
inefficiency on anisotropy is largest for the smallest i? max , and sta- 
tistically similar for the two larger values. Increasing the aperture 
to the virial radius or above increases the number of tracers near the 
virial radius where r2oo is estimated. Moreover, increasingly larger 
apertures will capture better the asymptotic value of the anisotropy 
profile A(r). In contrast, r p is less efficiently determined when 
fl max is increased from rsoo to rioo, while r„ has its worst ineffi- 
ciency for i? max = rtoo, with statistically equivalent values for the 
two smaller apertures. This might be due to the increasing fraction 
of unidentified interlopers, and/or to the presence of (sub) structures 
at larger radii. 

To assess the sensitivity of the MAMPOSSt technique to the 
number of tracers, besides the 500-particle selection, we have also 
considered samples of 100 particle tracers, randomly extracted 
from the same projections of the same 1 1 cosmological haloes. Re- 
sults of the MAMPOSSt analysis are listed at the bottom of Table 2 
and displayed in Fig. 3. Also in this case, for simplicity, we only 



© 2012 RAS, MNRAS 000, 000-000 



10 Mamon et al. 



o.io ; 
o cos ; 

0.00 [ 
-0.05 E 

-o.io : 

0.2 : 
0.1 

^ o.o : 
-o.i 

-0.2 : 
0.4 : 

o.2 ; 
^ o.o [ 

-0.2 F 
-0.4 



ft 




0.4 : 

o,2 ; 

0.0 | 
-0.2 : 
-0.4 : 



4f 



dHK KBM Clean (AT-100) 

Membership 



CO 

o 
cu 

% 1 



1.0 





i y 


r 200 














Oft/© 




cffF - 




o 


,<> 




V n 


1 




1 




true 



CO 

o 
a, 
s 
< 
s 



0.1 



o o, 
□ ^ 





o o 




°* 






° 




> r 




o 


o 



0.1 



true 




true 



1.0 




true 



Figure 4. Correlation of MAMPOSSt and true values of the 4 jointly-fit pa- 
rameters (Case Gen), with the T anisotropy profile, for each of the 3x11 
haloes with 500 tracers. Each panel corresponds to a different parameter, as 
labelled (units of radii are in h^ 1 Mpc). Different symbols identify differ- 
ent projections, ic-axis: black diamonds, y-axis: red squares, z-axis: blue 
circles. 



Figure 3. MAMPOSSt residuals, log(o/t), for the MAMPOSSt parame- 
ters (top: virial radius, 2nd panel: tracer scale radius, 3rd panel: DM scale 
radius, bottom: velocity anisotropy) for the different schemes of interloper 
removal (see text). The mean {dots) and dispersion (error bars) of log(o/i) 
are respectively illustrated as filled circles and error bars, for the 33 sam- 
ples of 500 ('Clean', dHK and KBM) and 100 ('Clean' (jV=100)) particles. 
Results for the anisotropy models Cst, T, and ML are shown left to right in 
green, blue, and red, respectively. 

display a limited set of results. When compared to the results for 
the 500-particle samples, there is no significant change in the av- 
erage values of the bias with which the different parameters are 
recovered, except for r p , where the bias becomes negative, while 
it was mostly positive for the 500-particle samples. The r p param- 
eter value underestimation is not very severe, however, ^ 25%. 
The efficiencies with which the different parameters are estimated 
are significantly affected by the reduction in number of tracers. The 
dispersion increases from ~ 10 to 15% for r2oo, from ~ 25 to 60% 
for TV, from ~ 60 to 100% for r p , and from ~ 20 to 33% for A 
and Aoo- There is no significant change in the dispersion for rp, 
but this was already extremely large for the 500-particle samples. 



3.4 Cases with constraints on parameters 

What is the effect of reducing the number of free parameters on the 
performance of the MAMPOSSt algorithm? To assess this point we 
consider several cases that reproduce what observers do in practice 
when faced with the problem of determining the internal dynamics 
of cosmological haloes. In all cases, we consider 500 particles se- 
lected within r2oo in each halo. We only apply the Clean interloper 



removal algorithm, and we only consider the NFW mass density 
model, for simplicity. 

Specifically, we consider the following Cases: 

A) General [Gen]: r2oo, r v , r p , and the anisotropy parameter 
(one among the following: A, rp, Aoo, depending on the anisotropy 
model considered) are all free MAMPOSSt parameters. This is the 
case considered so far. 

B) General with r v fitted outside MAMPOSSt [Split]: the free 
parameters are the same as in the Gen case, but r v is fitted out- 
side MAMPOSSt, via MLE. We thus split the minimization of the 
parameters into two parts. 

C) Known virial mass or radius [KVir]: r^oo is fixed and as- 
sumed to be exactly known, r p and the anisotropy parameter are 
free parameters in MAMPOSSt, r v is an external free parameter, 
as in the Split case. 

D) Estimated virial mass or radius [EVir]: similar to KVir, ex- 
cept that r2oo is not the true value, but the value estimated from 
the LOS aperture velocity dispersion (after interloper removal, see 
Appendix B). 

E) ACDM: r p is estimated from r2oo using the theoretical rela- 
tion between these two quantities provided by Maccio et al. (2008); 
r2oo and the anisotropy parameter are free parameters in MAM- 
POSSt, r„ is an external free parameter, as in the Split case. 

F) Mass follows Light [MfL]: r2oo and the anisotropy parameter 
are free parameters in MAMPOSSt, r v is an external free parame- 
ter, as in the Split case, and r p is assumed to be identical to r v . 

G) Tied Light and Mass [TLM]: r200 and the anisotropy param- 
eter are free parameters in MAMPOSSt, while r p and r v are as- 
sumed to be an identical free parameter. 

H) Isotropic [/3-iso]: A = 1 is assumed, r200,r p are free pa- 
rameters in MAMPOSSt, r v is an external free parameter, as in the 
Split case. 
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Figure 5. MAMPOSSt residuals as a function of true virial radius for the 4 
jointly-fit parameters (Case Gen), with the T anisotropy profile, for each 
of the 3x11 haloes with 500 tracers. Each panel corresponds to a different 
parameter, as labelled. Different symbols identify different projections, x- 
axis: black diamonds, y-axis: red squares, z-axis: blue circles. 



Table 3. Parameter inputs for the different cases of parameter constraints 



Case 




M(r) 






Tv 


r2oo 






Gen 


{R,Vz} 


{R,M 


{R, Vz} 


{R, v z } 


Split 


{R} 


{R, «4 


{R,Vz} 


{R, v z } 


KVir 


{R} 


known 


{R,vz} 


{R, Vz} 


EVir 


{R} 


from a z 


{R,Vz} 


{R, Vz} 


ACDM 


{R} 


{R, v z } 


from r-2oo 


{R, v z } 


MfL 


{R} 


{R, v,} 


Tv 


{R, Vz} 


TLM 


r P 


{R, v z } 


{R, v z } 


{R, Vz} 


/3-iso 


{R} 


{R, Vz} 


{R, vz} 





/3-MBM 


{R} 


{R, v z } 


{R, v z } 


rp=r P 


/3-HM 


{R} 


{R, v z } 


{R, vz} 


a + b d In u/d In r 



I) Anisotropy model a la MBM10 [/3-MBM]: ML anisotropy 
model with rp forced to be identical to r p ; r2oo,?> are free pa- 
rameters in MAMPOSSt, r v is an external free parameter, as in the 
Split case. 

J) Anisotropy linked to the mass density profile, using the 
anisotropy - slope relation of Hansen & Moore (2006) [/3-HM], 

Table 3 summarizes the different Cases, indicating, for each 



parameter, whether it is estimated from the full { R, v z } distribution 
with MAMPOSSt, the {R} distribution only, using standard MLE, 
or if it is fixed or linked to some other parameter. 

The Cases outlined above correspond to different observa- 
tional situations. Case Gen is the most general situation in which 
the observer ignores all the dynamical characteristics of the system, 
and it is the case analysed so far (see Table 2). We repeat it here for 
the sake of comparison with the other cases. 

Case Split is closer than case Gen to the typical observa- 
tional situation, since the number density profile of the tracers of 
the potential (stars, galaxies) is generally determined directly (e.g. 
by counting them in concentric radial annuli, if the system is as- 
sumed to be spherical). The advantage of Split relative to Gen is 
that any radial incompleteness that might affect the determination 
of r„ can be easily corrected for, if known, before running MAM- 
POSSt (see, e.g. Biviano & Poggianti 2009). However, as outlined 
in Sect 2.1, an interesting alternative is to include this incomplete- 
ness as in equation (20). Here, in all other cases (except TLM), r v 
has been determined outside MAMPOSSt. Anyway, the results (see 
below) obtained for the Gen and Split cases are very similar, hence 
it makes little difference if r„ is fitted within or outside MAM- 
POSSt, when the analysed sample (as in our case) does not suffer 
from radial-dependent incompleteness. 

In the KVir and EVir cases, r2oo is not a free parameter, but 
is fixed externally. In the KVir case, a perfect knowledge of r2oo 
is assumed, and in fact we take the true r2oo of the cosmological 
haloes. This case corresponds to the situation in which r2oo esti- 
mates are available from other data than the projected-phase space 
distribution of galaxies, e.g. from weak-lensing or X-ray observa- 
tions for clusters, although in the real world also these mass esti- 
mates will be affected by uncertainties, so KVir is a rather idealised 
case. Case EVir corresponds to the situation in which r2oo estimates 
are directly obtained from the Clean interloper removal scheme, as 
described in Appendix B. 

In the ACDM and MfL cases, r p is not a free parameter. In 
the ACDM case, we determine r p from r2oo using the relation of 
Maccio et al. (2008), which is based on the analysis of haloes in 
ACDM numerical simulations. Case MfL corresponds to the situa- 
tion in which the observer has good reasons (or a priori theoretical 
prejudice) to assume that 'Mass follows Light', i.e. that the tracer is 
spatially distributed like the mass. Therefore the observer first de- 
termines r„ from the distribution of the tracer and then makes the 



assumption r p 



In case TLM, the observer is unable to con- 



strain the tracer scale radius from its spatial distribution. This may 
happen when dealing with an incomplete sample with unknown 
incompleteness. Hence, r„ and r p are both determined from the 
dynamical analysis, assuming they are identical. 

In the /3-iso, /3-MBM, and /3-HM cases, the anisotropy is no 
longer a free parameter. The /3-iso case corresponds to the situa- 
tion where, for lack of better knowledge, the velocity anisotropy 
profile is assumed to be isotropic, A. = 1 (/? = 0) (like in, e.g. 
Biviano & Girardi 2003). Since the velocity anisotropy profiles of 
cluster-mass haloes have been shown to be well represented by an 
ML profile with rp — r p (MBM10, see also Fig. 1), it also makes 
sense to fix this anisotropy model in the fitting. This is the /3-MBM 
case. Finally, in the /3-HM case, we adopt an anisotropy profile that 
varies linearly with the logarithmic slope of the number density 
profile: 



/3(r) = a + b 



d In v 



(39) 



dlnr 

with a = —0.15 and b — —0.19, as was determined on a vari- 
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Figure 6. Same as Fig. 3 for the different cases of parameter constraints 
(see text and Table 3). As in Fig. 3, results for the anisotropy models Cst, 
T, and ML are shown, left to right, in green, blue, and red, respectively, 
and here the results for the /3-HM anisotropy model are shown in black. 
Open symbols refer to those parameters that are not free parameters in the 
MAMPOSSt analysis for the Case considered. 



ety of simulations by Hansen & Moore (2006). Note that Sparre 
& Hansen (2012) have argued that cosmological haloes are better 
described by this relation than by the attractor of Hansen, Juncher, 
& Sparre (2010), although Lemze et al. (2012) find that the validity 
of this relation is limited to the central regions of haloes. 

By fixing the anisotropy to determine the parameters of the 
mass profile, these cases are similar in spirit to mass inversion 
(Mamon & Boue 2010; Wolf et al. 2010), except that the former 
requires a parametric form for the mass profile, while the latter suf- 
fers from binning and extrapolation. 

So, while the anisotropy parameter is fixed in the /3-iso case, 
it changes according to r p (a free parameter here) in the /3-MBM 
case, and according to the logarithmic slope of the tracer number 
density profile (which is related to the free parameter r„) in the 
/3-HM case. 

In Table 4, we list the results of our analysis for the different 
observational cases for the 500-particle samples. These results are 
graphically displayed in Fig. 6. For the sake of comparison, we list 
and plot again here the results for the Gen case, already displayed 
in Table 2 and Fig. 3. 

The results for r2oo are essentially independent of the case 
considered. While the r2oo values for the EVir case are directly ob- 
tained from the Clean interloper removal scheme (App. B), they are 



measured with similar bias and inefficiency as for the Gen case. In 
other words, a z can be used to provide an estimate of r2oo with a 
comparable accuracy to that provided by MAMPOSSt. The advan- 
tage of MAMPOSSt is of course that it provides estimates for the 
other dynamical parameters at the same time. 

Fitting rv from the R distribution only, i.e. externally from 
MAMPOSSt, does not significantly alter the accuracy returned for 
this parameter. This means that incompleteness in the sample of 
tracers, if properly accounted for, is not a significant issue for 
MAMPOSSt. Moreover, this lifts our concern that the Split method 
does not find the same minimum for — In C as the Gen (joint) 
method (see end of Sect. 2.1). Using r p to predict r v (the TLM 
case) leads to larger bias and inefficiency on r v , but this occurs 
because the bias and inefficiency on r p are worse than those on r v . 

Also, the results for r p depend little on the case considered. 
The best results are obtained for the ACDM case, where r p is not a 
free parameter of the MAMPOSSt analysis, but it is estimated from 
r2oo using a theoretical relation. The good results obtained for the 
ACDM case are not surprising, given that r2oo is better constrained 
than r p in MAMPOSSt, and that the test haloes considered here are 
extracted from a ACDM cosmological simulation, and have there- 
fore similar properties to those used by Maccio et al. (2008) to es- 
tablish the mass-concentration (and hence the r2oo — r p ) relation. 

As for the anisotropy parameters, the ML model remains im- 
possible to constrain in all cases. For the Cst and T models, the 
bias and inefficiency do not depend strongly on the case considered. 
However, the o vs. t correlation values do depend on the case: a sig- 
nificant correlation is obtained for Aao for the Gen, Split, ACDM 
and TLM cases, but not for KVir, EVir, or MfL. For the ACDM and 
/3-MBM cases, a significant correlation is found even for the poorly 
constrained rp parameter, but note that in the /3-MBM case, rp is 
not a free parameter of the MAMPOSSt analysis, but is forced to 
match r p . 
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Table 4. MAMPOSSt results for 4 free parameters as well as several cases of constrained parameters 



Case 


0(r) 




T200 






r v 






r p 




anisotropy 








bias 


ineff. 


corr. 


bias 


ineff. 


COIT. 


bias 


ineff. 


corr. 


bias 


ineff. 


COIT. 


Gen 
Gen 
Gen 


Cst 
ML 
T 


0.004 

-0.003 
-0.006 


0.040 
0.040 
0.040 


0.909 
0.904 
0.903 


0.027 
0.024 
0.026 


0.102 
0.104 
0.103 


0.835 
0.832 
0.838 


0.032 
0.057 
0.039 


0.217 
0.229 
0.169 


0.578 
0.601 
0.709 


0.007 
-0.221 
0.007 


0.073 • 

0.887 

0.085 


-0.255 
-0.172 
0.621 


Split 
Split 
Split 


Cst 
ML 
T 


0.004 
-0.003 
-0.006 


0.040 
0.040 
0.040 


0.909 
0.899 
0.907 


0.032 
0.032 
0.032 


0.101 
0.101 
0.101 


0.833 
0.833 
0.833 


0.036 
0.044 
0.038 


0.216 
0.218 
0.173 


0.580 
0.591 
0.713 


0.010 

-0.072 
0.003 


0.074 
0.756 
0.083 


-0.277 
0.038 
0.627 


KVir 
KVir 
KVir 


Cst 
ML 
T 


E 


E 


E 


0.021 
0.021 
0.021 


0.103 
0.103 
0.103 


0.817 
0.817 
0.817 


0.017 
0.044 
0.040 


0.305 
0.244 
0.198 


0.266 
0.445 
0.638 


0.004 

-0.237 
-0.036 


0.084 
1.564 
0.143 


0.016 
0.221 
0.198 


EVir 
EVir 
EVir 


Cst 
ML 
T 


-0.000 
-0.000 
-0.000 


0.042 
0.042 
0.042 


0.884 
0.884 
0.884 


0.028 
0.028 
0.028 


0.103 
0.103 
0.103 


0.832 
0.832 
0.832 


0.041 
0.071 
0.028 


0.223 
0.201 
0.185 


0.596 
0.571 
0.714 


0.019 
-0.314 
-0.018 


0.064 ■ 

0.827 

0.086 


-0.155 
0.248 
0.431 


ACDM 
ACDM 
ACDM 


Cst 
ML 
T 


0.008 
0.002 
-0.004 


0.038 
0.039 
0.039 


0.906 
0.911 
0.903 


0.032 
0.032 
0.032 


0.101 
0.101 
0.101 


0.833 
0.833 
0.833 


0.070 
0.061 
0.061 


0.151 
0.154 
0.147 


0.674 
0.678 
0.675 


0.024 
-0.048 
0.040 


0.066 
0.797 
0.190 


0.470 
0.493 
0.497 


MfL 
MfL 
MfL 


Cst 
ML 
T 


0.005 
-0.001 
-0.003 


0.038 
0.039 
0.041 


0.901 
0.899 
0.887 


0.032 
0.032 
0.032 


0.101 
0.101 
0.101 


0.833 
0.833 
0.833 








0.002 
0.122 
0.011 


0.081 
0.906 
0.193 


0.127 
0.182 
0.194 


TLM 
TLM 
TLM 


Cst 
ML 
T 


0.003 
-0.006 
-0.003 


0.040 
0.046 
0.042 


0.908 
0.889 
0.906 








-0.035 
-0.009 
0.020 


0.276 
0.291 
0.213 


0.440 
0.422 
0.605 


-0.008 
0.116 

-0.024 


0.072 ■ 

0.715 

0.105 


-0.136 
0.183 
0.529 


/3-iso 

/3-MBM 

/3-HM 


Cst 
ML 

/3(P) 


0.003 
-0.005 
-0.003 


0.038 
0.040 
0.040 


0.914 
0.906 
0.899 


0.032 
0.032 
0.032 


0.101 
0.101 
0.101 


0.833 
0.833 
0.833 


-0.097 
0.024 
0.009 


0.215 
0.196 
0.222 


0.642 
0.651 
0.617 


-0.081 
0.044 


0.045 
0.455 


0.472 



Notes: These results are for 1 1 haloes of 500 particles, each observed along 3 axes out to the true value of r2oo in projection, with NFW density models and 
the 'Clean' interloper removal algorithm. Col. 1 : Case for MAMPOSSt analysis (see Table 3); col. 2: anisotropy model (Cst: = cst; ML: eq. [37], Mamon & 
Lokas 2005; T: eq. [38], adapted from Tiret et al. 2007); cols. 3-5: virial radius; cols. 6-8: tracer scale radius; cols. 9—11: dark matter scale radius; cols. 12-14: 
velocity anisotropy (i.e., A for the Cst model, rg for the ML model, and *4oo for the T model). The columns 'bias' and 'ineff.' respectively provide the mean 
and standard deviation (both computed with the biweight technique) of log(o/i), while columns 'corr.' list the Spearman rank correlation coefficients between 
the true values and MAMPOSSt-recovered ones (values in boldface indicate significant correlations between o and t values at the ^ 0.99 confidence level). 
Values in italics indicate parameters that are not free in the MAMPOSSt analysis for the Case considered. 



It is surprising that MAMPOSSt obtains worse inefficiencies 
in the KVir case where T2oo is assumed perfectly known in compar- 
ison to the Gen, Split or EVir cases. The reason for this is probably 
related to the halo triaxialities. When we observe a halo along a 
given line-of-sight, we are sampling its velocity distribution only 
along that line-of-sight. If the components of the halo velocity dis- 
tribution along different axes are different, this difference will be 
reflected in our results. 

This point is demonstrated in Fig. 7. There, we show the cor- 
relation between the ratio of r 2 oo measured by MAMPOSSt to the 
true value measured in 3D versus the ratio of velocity dispersions 
within the true virial sphere, measured along the projection axis 
over that measured globally. The very high correlation (Spearman 
rank: r = 0.82) shows that the error on r2oo is related to the veloc- 
ity dispersion ratio above, which directly measures the triaxiality 
of the halo (without mixing with the effects of interlopers beyond 
the virial sphere). 



3.5 Stacked haloes 

To extract all the possible information from the available data, it is a 
common practice to construct stacked samples (e.g. Carlberg et al. 



1997; Biviano & Girardi 2003; Katgert et al. 2004; Biviano & Pog- 
gianti 2009). In the case of clusters of galaxies, this is done by scal- 
ing the galaxy distances from their cluster centres and the galaxy 
velocities with respect to their cluster mean velocity, by their clus- 
ter r 2 oo and v 2 <m, respectively, where v 2 <m = (GM200A200) 1 ' 2 ■ 
To do this, prior knowledge of the individual cluster r2oo values 
is needed. Trying to mimic the observational situation as close as 
possible, we stack our haloes using the values of r 2 oo directly ob- 
tained from the Clean interloper removal scheme, as described in 
Appendix B, i.e. the values used in what we called the EVir case 
in Sect. 3.4. We thus build three stacked haloes, one for each pro- 
jection axis, from the eleven 500-particle haloes, each first passed 
through the Clean interloper removal scheme (app. B). 

The stacks thus created from the 500-particle (respec- 
tively 100-particle) samples contain 5248, 5260, 5213 (respec- 
tively 1061, 1065, 1048) DM particles along the a>, y-, and z-axis, 
respectively. 

Since we need to fix the r2oo values of the haloes before stack- 
ing, r2oo cannot be a free parameter of the MAMPOSSt analysis. 
Moreover, in an effort to mimick the observational procedure, we 
estimate the individual r v values of the 1 1 haloes also before stack- 
ing, and then we take the biweight average of these values as rep- 
resentative of the rv of the stack. In the real world, this is done 
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Table 5. MAMPOSSt results for the stacked haloes 



Projection 


0{r) 


»~200 


TV 


r 


p 


anisotropy 










stack 


mean 


stack 


mean 


11 x 500 


X 


Cst 


1.104 


0.208 


0.246 


0.227 


1.195 


1.230 


y 


Cst 


1.050 


0.186 


0.199 


0.292 


1.184 


1.253 


z 


Cst 


1.018 


0.136 


0.146 


0.260 


1.130 


1.240 


bias 




-0.014 


-0.166 


-0.129 


-0.003 


-0.008 


0.011 


inefficiency 




0.021 


0.109 


0.134 


0.065 


0.014 


0.005 


X 


ML 


1.104 


0.208 


0.219 


0.328 


0.801 


0.045 


y 


ML 


1.050 


0.186 


0.203 


0.306 


0.429 


0.049 


z 


ML 


1.018 


0.136 


0.154 


0.259 


0.494 


0.256 


bias 




-0.014 


-0.166 


-0.120 


0.060 


0.305 


-0. 743 


inefficiency 




0.021 


0.109 


0.090 


0.061 


0.160 


0.446 


X 


T 


1.104 


0.208 


0.283 


0.337 


1.565 


1.330 


!) 


T 


1.050 


0.186 


0.235 


0.256 


1.509 


1.330 


Z 


T 


1.018 


0.136 


0.180 


0.252 


1.459 


1.427 


bias 




-0.014 


-0.166 


-0.054 


-0.010 


0.018 


-0.038 


inefficiency 




0.021 


0.109 


0.116 


0.075 


0.018 


0.018 


11 X 100 


X 


Cst 


1.107 


0.150 


0.143 


0.160 


1.101 


1.164 


y 


Cst 


1.016 


0.190 


0.144 


0.171 


1.224 


1.402 


z 


Cst 


1.038 


0.142 


0.146 


0.210 


1.151 


1.181 


bias 




-0.017 


-0.238 


-0.256 


-0.169 


-0.019 


-0.014 


inefficiency 




0.022 


0.075 


0.050 


0.070 


0.027 


0.048 


X 


ML 


1.107 


0.150 


0.155 


0.178 


0.515 


0.118 


y 


ML 


1.016 


0.190 


0.168 


0.212 


0.093 


0.051 


z 


ML 


1.038 


0.142 


0.163 


0.228 


0.263 


0.226 


bias 




-0.017 


-0.238 


-0.205 


-0.100 


-0.043 


-0.368 


inefficiency 




0.022 


0.075 


0.021 


0.064 


0.439 


0.382 


X 


T 


1.107 


0.150 


0.189 


0.183 


1.515 


1.506 


.'/ 


T 


1.016 


0.190 


0.164 


0.202 


1.464 


1.424 


z 


T 


1.038 


0.142 


0.183 


1.424 


1.554 


1.276 


bias 




-0.017 


-0.238 


-0.158 


-0.130 


0.018 


-0.014 


inefficiency 




0.022 


0.075 


0.036 


0.025 


0.015 


0.042 



Notes: Col 1: viewing axis; col. 2: anisotropy model; col. 3: virial radius (bi- 
weight mean over 1 1 haloes from Clean interloper removal scheme); col. 4: 
tracer scale radius (biweight mean of MAMPOSSt EVir estimate for the 
11 haloes); col. 5: DM scale radius (MAMPOSSt); col. 6: same (biweight 
mean of MAMPOSSt EVir estimate for the 1 1 haloes); col. 7: anisotropy 
parameter (A for Cst, rp for ML and Aaa for T) from MAMPOSSt; col. 8: 
same (biweight mean of MAMPOSSt EVir estimate for the 1 1 haloes). Val- 
ues that are not directly obtained through MAMPOSSt on the stacked halo 
but from the mean of the individual haloes are shown in italics. The bi- 
ases and inefficiencies are respectively computed from biweight means and 
gapper standard deviations of the 3 values (see Beers et al. 1990). 



because individual haloes often suffer from different spectroscopic 
incompleteness levels. More precisely, we take the r u estimates ob- 
tained in the EVir case (see Sect. 3.4). The remaining two parame- 
ters, r p and the anisotropy parameter are estimated via the MAM- 
POSSt analysis. Results for the six stacked haloes are listed in Ta- 
ble 5 and displayed in Fig. 8. We use the gapper estimate of dis- 
persion (Wainer & Thissen 1976, see Beers et al. 1990) given our 
small sample (N = 3). 

In Fig. 9, we show 1-, 2-, and 3-cr contours in the logr p — 
log Aoo plane (for the T anisotropy model). In this figure, we indi- 
cate with a cross the expected solution for the stacked haloes. This 
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Figure 8. MAMPOSSt residuals for stacked haloes (labelled 'S500' and 
'S100', when built from the 11 haloes, respectively sampled with 500 and 
100 particles), and biweight mean of the values obtained for the 1 1 haloes 
in the EVir case (see Sect. 3.4, labelled 'M500' and 'Ml 00', for the 500 
and 100 particle haloes, respectively). Different symbols identify different 
projections, rr-axis: diamonds, j/-axis: squares, z-axis: circles. In the upper- 
left panel (r p ), the green, blue, and red symbols are the results obtained 
using, respectively, the Cst, T, and ML anisotropy models. 



is the biweight-average (with its error) of the true r p and Aoo val- 
ues of the 1 1 individual haloes from which the stacked haloes are 
constructed. These values are also listed in the last line of Table 1. 
The average values of r_2 and rp are identical, in substantial agree- 
ment with what was found by MBM10. Note also that the average 
value of Aoo corresponds to j3oo = 0.5, meaning that the average 
ML and T anisotropy models are the same. 

Fig. 9 shows that the expected solution for the stack sample 
is always within the 1-ct contour of the MAMPOSSt result for the 
T model. While not shown, this is also the case for the Cst model, 
but not for the ML model, where rp remains essentially uncon- 
strained, as was the case for the individual haloes. The estimates of 
the parameter r p do not depend on the assumed anisotropy model. 
They are biased low for the stacks built from the 100-particle sam- 
ples, even if not significantly so (see Fig. 9). Interestingly, while the 
bias is stronger for the 100-particle stack than for the 500-particle 
one, the dispersion of the recovered values appears to be lower. The 
dispersion is instead higher, as expected, for the anisotropy param- 
eters, except for Aoo where the dispersion is comparable for the 
500-particle and the 100-particle stacks. However, with only three 
results per anisotropy model for each Stack sample, these differ- 
ences in bias and dispersions do not appear to be statistically sig- 
nificant. 

In Table 5, we also list the biweight average values of the 1 1 
individual haloes along each projection, obtained in the EVir case 
(see Sect. 3.4), and these values are also plotted in Fig. 8 (labelled 
'M') together with the best-fit solutions for the stacked haloes (la- 
belled 'S'). Comparing the best-fit values of r p and the anisotropy 
parameters obtained from the stacked haloes with those found from 
the mean of the individual 1 1 haloes (along the same projection), 
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Figure 9. MAMPOSSt confidence contours in the log r p — log Aoo plane 
for the stacked haloes constructed from single-axis projections of 1 1 haloes 
from the 500-particle (left panels) and 100-particle (right panels) samples. 
The x-axis, j/-axis, and 2-axis stacks are shown in the top, middle, and bot- 
tom panels, respectively. The contours are drawn at 1-, 2- and 3-cr. The 
best-fit solutions are indicated by the filled circles, while the crosses rep- 
resent the biweight average of the true values of the parameters for the 1 1 
haloes. The cross lengths across the x- and j/-axes represent the errors on 
the two parameter averages (see Table 1). 



it is difficult to decide whether it is better to run MAMPOSSt on 
the stacked halo or to adopt means of the individual MAMPOSSt 
estimates. A close inspection of Table 5 and Fig, 8 indicates that, 
for iV = 500, stacks are more accurate for estimating rp and Aoo, 
while means of individual MAMPOSSt estimates appear more ac- 
curate for estimating r p and A. But the statistics are poor. In the 
real observational situation, it therefore seems advisable to con- 
sider both the average of the results of the individual haloes and the 
result obtained for the stacked sample. 

3.6 A large halo 

It is interesting to consider whether we can achieve the same accu- 
racy using a single halo, but with a total number of tracer particles 
similar to that of the stack. We extract ~ 5000 particles within r2oo 
for the more massive halo in our sample, 5726, remove interlopers 
with the Clean method, and run MAMPOSSt for the general case, 
adopting the NFW mass density and different anisotropy models. 

Table 6 compares the results for the 3 cones of 500 particles 
(counted before interloper removal) and the 3 with ~5000 parti- 
cles. We use again the gapper dispersion, given our small sample 
(N = 3). As expected, the biases and inefficiencies are gener- 
ally reduced when selecting all (~5000) particles instead of ran- 
domly selecting 500. For the T anisotropy model, the inefficiencies 
in log(o/t) are reduced by 0.1 dex for r p and anisotropy, 0.15 dex 
for r v , but virtually unchanged for r2oo- Surprisingly, for Cst and 
ML anisotropy, r ^oo is reached with worse inefficiency with the full 
5000-particle sample in comparison with the 500-particle one. 

However, one should not over-interpret these comparisons. In- 
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Table 6. MAMPOSSt results for the most massive halo, 5726 



Projection 


JV 


/J(r) 




Tv 


rp 


anisotropy 


X 


4529 


Cst 


1.538 


0.499 


0.262 


1.053 


II 


4520 


Cst 


1.747 


0.369 


0.352 


1.122 


z 


4662 


Cst 


1.642 


0.318 


0.197 


1.252 


bias 






-0.005 


-0.024 


-0.190 


-0.096 


inefficiency 






0.033 


0.116 


0.149 


0.044 


X 


4529 


ML 


1.505 


0.499 


0.308 


0.568 


11 


4520 


ML 


1.711 


0.370 


0.352 


0.638 


z 


4662 


ML 


1.621 


0.318 


0.190 


0.301 


bias 






-0.013 


-0.023 


-0.142 


1.080 


inefficiency 






0.033 


0.115 


0.158 


0.192 


X 


4529 


T 


1.454 


0.502 


0.447 


2.289 


V 


4520 


T 


1.700 


0.371 


0.425 


1.666 


z 


4662 


T 


1.623 


0.320 


0.236 


1.705 


bias 






-0.017 


-0.021 


0.030 


-0.116 


inefficiency 






0.040 


0.115 


0.164 


0.082 


X 


489 


Cst 


1.616 


0.492 


0.336 


0.887 


y 


460 


Cst 


1.773 


0.398 


0.340 


1.156 


z 


476 


Cst 


1.688 


0.260 


0.187 


1.188 


bias 






0.008 


-0.036 


-0.081 


-0.084 


inefficiency 






0.024 


0.164 


0.154 


0.075 


X 


489 


ML 


1.554 


0.494 


0.483 


0.764 


V 


460 


ML 


1.739 


0.399 


0.335 


0.562 


z 


476 


ML 


1.672 


0.261 


0.195 


0.310 


bias 






-0.001 


-0.034 


-0.107 


1.015 


inefficiency 






0.029 


0.164 


0.232 


0.231 


X 


489 


T 


1.454 


0.496 


0.676 


3.454 


V 


460 


T 


1.724 


0.401 


0.412 


1.735 


z 


476 


T 


1.677 


0.261 


0.247 


1.829 


bias 






0.010 


-0.033 


0.003 


-0.092 


inefficiency 






0.044 


0.165 


0.258 


0.177 



Notes: The results are for Case Gen. Col. 1: viewing axis; col. 2: number of 
particles after interloper removal with the Clean method; col. 3: anisotropy 
model; col. 4: virial radius; col. 5: tracer scale radius; col. 6: DM scale 
radius; col. 7: anisotropy parameter (A for Cst, rp for ML and Aoo for 
T); The biases and inefficiencies are respectively computed from biweight 
means and gapper standard deviations of the 3 values (following Beers et al. 
1990). 

deed, our Monte-Carlo tests on 10 000 random samples of 3 objects 
arising from a Gaussian distribution indicate that while the gapper 
standard deviation is unbiased even for as few as 3 objects, it has 
an inefficiency of 0.53 times its true value for N — 3. In other 
words, an inefficiency of log(o/i) of 0.1 would have roughly 0.05 
accuracy. Therefore, any improvement or worsening of inefficiency 
in log(o/t) of less than a factor of 2 is clearly not statistically sig- 
nificant. Nevertheless, the fact that most inefficiencies are reduced 
when increasing the sample size from 500 to 5000 is much more 
significant, and is indeed what is qualitatively expected. However, 
the expected quantitative improvement of log VT0 = 0.5 does not 
appear to be attained for the inefficiency of log(o/t). 

The bias and inefficiencies for the 5000-particle sample are 
typically higher for the single halo than for the 11 x 500 particle 
stack (with just 10% more particles, although the extraction there 
was for case EVir instead of Gen). But again, the inefficiencies are 
within a factor 2, so these changes are individually not statistically 
significant. But the combination of all of them appears to be statis- 
tically significant. 
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The worse inefficiencies for the large halo in comparison with 
the stacked halo of similar number of tracers can be due to either 
our assumption of Gaussian velocity distributions (see Sect. 2.2) or 
to the triaxiality of the halo velocity ellipsoids. The velocity dis- 
tribution of the stack halo approaches a Gaussian by the central 
limit theorem, so our Gaussian assumption is better verified in a 
stack sample than in individual haloes. As a matter of fact, a num- 
ber of studies have shown the 3D velocity distributions of individ- 
ual haloes to deviate from Gaussianity (Wojtak et al. 2005; Hansen 
et al. 2006). But perhaps the dominant effect is that of triaxiality. 
As discussed at the end of Sect. 3.4, the results obtained by MAM- 
POSSt for r2oo are influenced by the choice of the projection axis, 
and in the present case, the r2oo values obtained for 5726 along the 
x, y, and z projections when using 5000 particles, are ordered in 
the same way as those obtained when using 500 particles. 



4 CONCLUSIONS AND DISCUSSION 

We have presented the formalism for a new method, called MAM- 
POSSt, for the determination of the mass and anisotropy profiles of 
spherical systems, assumed to be in dynamical equilibrium, in par- 
ticular, as stationary systems with no streaming motions. In MAM- 
POSSt, the distribution of tracers in projected phase space are fit 
by the predicted distribution arising from assumed 3D radial pro- 
files of the tracer density, total mass, and velocity anisotropy, as 
well from an assumed family of shapes for the radial and tangential 
components of the 3D velocity distribution. 

MAMPOSSt has several important advantages over other 
methods for mass/anisotropy modeling: 

1) MAMPOSSt does not involve binning the data, hence its re- 
sults are independent of the choice of such radial bins, in contrast 
with methods based upon velocity moments, including anisotropy 
and mass inversions; 

2) MAMPOSSt does not involve interpolations and extrapola- 
tions of binned radial profiles of the observables, again in contrast 
with anisotropy and mass inversions; 

3) MAMPOSSt does not require differentiation of the ob- 
served projected pressure, £ (-R) a'i{K), once more in contrast with 
anisotropy and mass inversions; 

4) MAMPOSSt extracts accurate constraints on the velocity 
anisotropy, in contrast with methods that assume Gaussian LOS 
velocity distributions; 

5) MAMPOSSt is very fast, as it involves a single integral (for 
popular f3(r) profiles). Indeed, the calculations for the best fit pa- 
rameters and their marginal and correlated distributions through 
MCMC for a 500-tracer system (as displayed in Fig. 2) require 
roughly 4 minutes of CPU time on a standard desktop personal 
computer. By contrast, distribution function methods involve triple 
integrals, and therefore take typically 1000 times longer to run, i.e. 
a few days for 500-tracer systems. Orbit modeling techniques are 
even slower and cannot properly probe parameter space. 

In this work, we have extensively tested MAMPOSSt, for the 
case of Gaussian 3D velocity distributions, on a set of 1 1 cluster- 
mass haloes extracted from cosmological simulations. The results 
of these tests indicate that, for systems with 500 velocities, MAM- 
POSSt provides essentially unbiased estimates of the relevant mass 
and velocity anisotropy profile parameters, with inefficiencies of 
10% for the virial radius, 20% for the constant or outer velocity 
anisotropy, 27% for the tracer scale radius, but as high as 48% 




0.8 0.9 1 1.1 1.2 



Figure 10. Comparison of virial radius errors from two schemes: MAM- 
POSSt vs. Clean interloper removal scheme (appendix B), for the 3x11 
haloes. The red line corresponds to equality. 



for the scale radius of the DM. However, MAMPOSSt seems un- 
able to set constraints on the radius of transition of a gently rising 
anisotropy model such as ML. 

We have noted that the results of MAMPOSSt are similar 
when we inserted the virial radius among the parameters to be 
jointly fit or when we derived the virial radius using the new in- 
terloper removal scheme that we presented in Appendix B. We 
found that this interloper-removal algorithm performs as well or 
better than other methods. As shown in Fig. 10, the errors on r2oo 
from MAMPOSSt are highly correlated with those obtained from 
the Clean procedure. Hence, the new interloper removal scheme 
produces extremely fast estimates of the virial radius, based on the 
LOS velocity dispersion, that are very close to those obtained with 
the full MAMPOSSt procedure. 

The correlation highlighted in Fig. 10 suggests that the qual- 
ity of MAMPOSSt estimates is limited by the accuracy of the in- 
terloper removal. However, we also saw (Fig. 7) that the error on 
log r2oo is correlated with the ratio of velocity dispersions within 
the true virial sphere, measured along the projection axis over that 
measured globally. This second correlation indicates that our capa- 
bility of recovering the true values of r2oo is limited by the effects 
of triaxiality, i.e. by the fact that we only observe one component 
of the 3-dimensional velocity distribution. The effect of triaxiality 
may be even more important than the uncertainties in the interloper 
cleaning process. 

Our additional tests lead to the following conclusions: 

1) Best results are obtained setting all 4 parameters free and us- 
ing the T anisotropy model. 

2) Setting the virial radius to the true value (e.g. measured 
through other techniques, such as X-rays or lensing) leads to worse 
results on the other parameters than when the virial radius is an 
additional free parameter. We argue that this is caused by the triax- 
iality of the haloes (Fig. 7). 

3) MAMPOSSt works surprisingly well with samples of 100 
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tracers instead of 500, only somewhat better with samples of 5000 
tracers, even better with stacked haloes of 11 x 500 tracers, simi- 
lar to the mean of these individual tracers, but not well as expected 
from N" 1 / 2 arguments. 

4) MAMPOSSt is not very sensitive to the chosen aperture, 
when this is to within «35% of the true value of r2oo- 

Table 7 compares the accuracy of MAMPOSSt with that of 
the dispersion-kurtosis method of Lokas (2002, tested by Sanchis, 
Lokas, & Mamon 2004) and the DF method of Wojtak et al. (2009). 
In comparison with the dispersion-kurtosis method, MAMPOSSt 
does slightly less well on the inefficiencies of the virial radius and 
worse on that of the DM scale radius (but better on the biases 
of both radii), but MAMPOSSt performs better for the constant 
anisotropy (both bias and inefficiency). Moreover, MAMPOSSt has 
the advantage of being flexible enough to infer the outer anisotropy, 
while the dispersion-kurtosis method of Lokas (2002) is limited to 
constant anisotropy (but this limitation to the dispersion-kurtosis 
method has been recently lifted by Richardson & Fairbairn 2012). 
Comparing with the much slower DF method, MAMPOSSt does 
somewhat worse on the inefficiency on the scale radius, with the 
same bias, yet MAMPOSSt is respectively slightly (not signifi- 
cantly) and much better on the inefficiencies on the virial radius 
and outer anisotropy, with less bias on the former. Note, however, 
that we forced here the inner velocities to be isotropic, while these 
were free parameters in the DF method and that the DF method 
computes the outer anisotropies at 5 r s instead of the asymptotic 
value. Nevertheless, we found (Table 4) that MAMPOSSt always 
performs as well or better when all parameters are set free, and we 
would expect the same for the DF method. Moreover, we see no 
reason why the inefficiency of the anisotropy measured near the 
virial radius should be worse than that of its asymptotic value. 

Despite the very good results of our tests, it must be noted that 
although the Gaussian MAMPOSSt recovers the 2nd velocity mo- 
ment (eq. [30]) , it does not recover the fourth moment. Indeed, tak- 
ing the fourth moment of the velocity distribution of equation (28) 
leads to 



X(R)vj(R) = 



v\ g(R, v z ) dv z 



v{r) 



r dr 



n J R a z (R,r) Vr 2 - R 2 



2a*(R,r) 



= 6 



p2 

1 



dv z 
r dr 



Vr 2 - R 2 



(40) 



where the second equality of equation (40) is obtained after re- 
versing the order of integration, while the final equality uses equa- 
tion (26). Equation (40) then gives the LOS velocity kurtosis excess 



Kz(R) 



vl(R) 
oi{R) 



(41) 



Unfortunately, the expression for the fourth velocity moment 
in equation (40) differs from the expression found by Lokas (2002) 
from the 4th order Jeans equation in the case of f3 = est. At large 
projected radius (R ~ r v ), k z estimated by equations (40) and (41) 
turns out to be nearly independent of (3. This shows the limit of the 
Gaussian approximation for the 3D velocity distribution. 

Nevertheless, despite the poor adequacy of the Gaussian ap- 
proximation, which does not produce the correct LOS kurtosis pro- 



files, MAMPOSSt with Gaussian velocities performs quite well 
according to our numerous tests. We are therefore confident that 
the method is mature for being used on real data-sets. And we are 
preparing analyses on several scales: dwarf spheroidal galaxies, gi- 
ant ellipticals (traced by planetary nebulae), groups and clusters of 
galaxies. We also plan to test MAMPOSSt on elliptical galaxies 
formed in dissipative mergers. 

One should note that in ACDM haloes, the 3D velocity distri- 
bution, measured in shells, is not Gaussian (Wojtak et al. 2005), but 
resembles the q-Gaussian distribution, which generalizes the Gaus- 
sian distribution in the same way as Tsallis (1988) entropy gen- 
eralizes Boltzmann-Gibbs entropy. The q-Gaussian (often called 
Tsallis) fits well a host of velocity distributions found in single- 
component dissipationless self-gravitating systems, with an index 
that varies linearly with the slope of the density profile (Hansen 
et al. 2006). We are thus preparing a generalization of MAMPOSSt 
to the Tsallis distribution of 3D velocities (appropriately combining 
radial and tangential terms), and expect this to perform even better 
than the MAMPOSSt-Gaussian algorithm. 
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Table 7. Comparison of the biases and inefficiencies of MAMPOSSt with those of other methods 



Method 


Reference 


Particles 










A 








bias 


ineff. 


bias 


ineff. 


bias 


ineff. 


bias 


ineff. 


Dispersion-kurtosis 


Sanchis et al. (2004) 


400 


-0.023 


0.033 


0.080 


0.240 


-0.040 


0.110 






MAMPOSSt /3=cst 


this article 


400 


0.004 


0.040 


-0.030 


0.289 


-0.010 


0.069 






Distribution function 


Wojtak et al. (2009) 


300 


-0.026 


0.048 


0.017 


0.169 






-0.002 


0.212 


MAMPOSSt /3-T 


this article 


300 


-0.003 


0.044 


0.017 


0.198 






-0.016 


0.116 



Notes: The biases and inefficiencies on virial masses for the dispersion-kurtosis and distribution function methods were divided by 3 to convert to biases and 
inefficiencies for virial radii. All analyses are for case TLM, and for MAMPOSSt the Clean interloper removal scheme is used. Dispersion-kurtosis log(o/t) 
for the DM scale radius were taken from the analogous values that Sanchis et al. (2004) obtained for the concentration parameter (basically assuming that 
in comparison, the errors on the virial radius are small). For the distribution function method, values of log(o/i) for the outer anisotropy were measured at 
roughly the virial radius, while the MAMPOSSt comparison (last row) is for the T anisotropy model (eq. [38]). The biases and inefficiencies are estimated 
with the biweight mean and dispersion, except for the dispersion-kurtosis method, where they are taken from the 2nd line of Table 3 of Sanchis et al.. Values 
where MAMPOSSt is better are highlighted in bold. 
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APPENDIX A: KERNELS FOR DETERMINING THE RADIAL VELOCITY DISPERSION PROFILE FROM THE JEANS 
EQUATION 

The radial squared velocity dispersion profile of equation (9) can be written (van der Marel 1994; Mamon & Lokas 2005) 



er r (r) = 



K(r)u(r) J r 



K(s) v(s) ^-t- as 



where K(r) is the solution to d In K/d In r — 2 /3(r), i.e. 



KM 



exp 



dt 



(Al) 



(A2) 



For simple anisotropy models (not all are used in this work), 



2 . 2 

r +r 8 



P = est, 
P = Pom 



r 2 + ri 



r + : 



P = Pml = - 



2 r + r fi ' 



(r + r p ) 



2ri, 



P = Pt = P a 



r + rp ' 



K{r) = ^ 



r + rg' 



(A3) 



exp 



1/3 



a + b 



7o r- + 7oo r 
r + r~ 



P = Pums (r < rp), 

P = Pbms (r ^ rp), 
d In v 



P = a + b 



2a / i t l/m\ o 

r exp ( — 4 bmr ' J p = 



dlnr 

d In 


dlnr 



(Hansen & Moore / Zhao) 



(Hansen & Moore / Einasto) 



where pou is the Osipkov-Merritt anisotropy model (Osipkov 1979; Merritt 1985), /3ml is the Mamon-Lokas anisotropy model (Mamon & 
Lokas 2005), /3|f n is a generalization of the Mamon-Lokas model by Tiret et al. (2007), /3dms = Min [l, {r/rp) 1 ^^ is another model that 
fits well the anisotropy of ACDM haloes with rp — r v (Diemand, Moore, & Stadel 2004), while the last two rows describe the anisotropy 
that varies linearly with the slope of the density profile (Hansen & Moore 2006), since the kernel depends on the precise form of the density 
profile. The first of these last two rows gives the kernel for Zhao's (1996) general density profile 



p(r) oc 



where 70 and 700 are the inner and outer logarithmic slopes of the density profile (i.e. 70 = — 1 and j a 



(A4) 



-3 for the NFW model), and 



-[(2 + 700)/ (2 + 70)] r_2 is the radius of slope 7 = (70 + 700 )/2, while r„2 is the radius of slope —2. The last row is for the Einasto 



(1965) density profile 



p(r) oc exp 



-2 m 



r 
r-2 



l/m 



(A5) 



which, for m ~ 6, is an excellent fit to the density profiles of simulated ACDM haloes (Navarro et al. 2004), while m = 5 is an excellent fit 
to the dark matter density profiles of haloes in hydrodynamical ACDM cosmological simulations (MBM10). 



APPENDIX B: INTERLOPER REMOVAL 

This appendix describes how we remove obvious (high absolute line-of-sight velocity) interlopers. 

We analyze the DM particles in the z — output of the hydrodynamical cosmological simulation of Borgani et al. (2004), to which 
we added the Hubble flow, placed an observer at 90 h^ 1 Mpc from each of the 105 most massive haloes (including the most irregular ones), 
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in each of the three cartesian directions. We limited these 3 x 105 = 315 haloes to projected radii within i? max = X max r 2 oo from the 
bary centre of the Friends-of-Friends halo identified in real space, where Jf max = 0.7, 1 or 1.4, and within line-of-sight velocities within 4 
times the true circular velocity at the V2oo from the true velocity centre (see MBM10 for details). The choice for the three values of X mn is 
to consider cases where the observers guess incorrectly the virial radius when they select their galaxies. We only used the 274 haloes out of 
315 that had at least 500 particles with R < r 2 oo- 6 

We chose TV particles at random in the projected phase space, so that we ended up with approximately n — 500 particles with R < r 2 oo, 
i.e. N = 500 for X max = 1 or N = 500 M™ FW (X m ^r 2a0 )/M™ FW (r 2a0 ) = 386 and 620 for X max = 0.7 and 1.4, respectively, where 
we assumed a concentration c = r 2 oo /f-2 = 4, typical of the haloes of our simulation (MBM10). The projected mass of the NFW model is 
M P (R) = M(r_ 2 ) M p (l£/r_ 2 ), where r_2 is the radius where the logarithmic of the density profile is —2 and the dimensionless projected 
mass is (Lokas & Mamon 2001, first derived by Bartelmann 1996 in a slightly longer and more computer intensive form) 



cosh _1 (l/X)/Vl -X 2 +ki(X/2) toiX < 1 , 

M P (X) = , n \ /n { l-ln2 forX = l, (Bl) 

:os- 1 (l/X)/VX 2 - 1 + ln(X/2) forX > 1 . 



In 2 - 1/2 



Our algorithm for interloper rejection, which is Bayesian as it assumes what we know about ACDM haloes, namely that they approxi- 
mately follow NFW density profiles with increasingly radial orbits, goes as follows. 

i) On first pass, we apply the gapper (or weighted gap) technique in order to split multimodal distributions of the distribution of velocities 
Vi, and identify the peak in velocity space closest to the input mean velocity of the halo, which we assume to know. After sorting the velocity 
offsets, we compute weighted gaps 

Qi = [i(n-i) (vi+i -Vi)] 1 ' 2 (B2) 

for 1 ^ i < n. We then check if the largest value of the dimensionless gap C«/MidMean(<5) is greater than some threshold, commonly 
called C (MidMean is the arithmetic mean within the quartiles of the distribution). Among all the subsamples separated by a dimensionless 
gap ^ C, we keep only that one closest to the input mean velocity of the halo. The gapper technique (Wainer & Thissen 1976) was first used 
in astronomy in the ROSTAT package 7 , which recommends C = 2.25 and was first applied to detect multimodal populations in clusters of 
galaxies by Girardi et al. (1993), who used C = 4. This is also the value we use here. 

ii) Using the velocities of the particles in the subsample identified with the gapper technique, we compute, in a first step, the global 
velocity dispersion of our selected particles, using the robust median absolute deviation (MAD), with <t v ,mad — MAD(i>)/0.6745 where 
MAD = Median|t> - Median(u)] (e.g. Beers et al. 1990). 

iii) We then estimate the virial velocity , v^ st = « c i rc (r^ st ) from the relation 

(rv)V _ [l/Aipfrv)] / rv 27vRE(R)a 2 z (R)dR 
GM(r v )/r v 

rr v poo , 

) I RdR L K >(i>i)« r)M{r) T (B3) 

XdX f K,(%, ^) p(x) M(x) ^ , (B4) 



M{r v )M p (r vi 



M(c)M p (c) J 



where p{x) = /o(r_ 2 x) / [M(r_ 2 )/(47rr?. 2 )], M(x) = M(r_ 2 a;)/M(r- 2 ), and M p (x) = M p (r_ 2 x)/M(r- 2 ) are the dimensionless 
radial profiles of density, mass and projected mass, respectively. Equation (B3) was derived using the relation (Mamon & Lokas 2005) 

X(R)a 2 z (R) = 2G J°°K fi ^) p(r)M(r) ^ (B5) 

with Kp a kernel that has been derived by Mamon & Lokas for simple anisotropy profiles. In equation (B4), we assume ML anisotropy 
(eq. [37]) with — r_2 as found in a cosmological simulation by MBM10. We then derive v v = v c i TC (r v ) from equation (B4). We adopt 
the NFW model, with dimensionless mass density p(x) = x^ 1 {x + l)~ 2 /(ln 2 — 1/2), dimensionless mass M(x) — [\n(x + 1) — x/(x + 
l)]/(ln 2 — 1/2) and dimensionless projected mass density M P (X) given in equation (Bl). With the Mamon & Lokas anisotropy profile, 
and rp = r_ 2 , the aperture velocity dispersion is well approximated (to better than 0.1% relative accuracy in the range 1 ^ c < 32) by 



3 



^a»(log c) 



(B6) 



where ao = —0.1197, ai = —0.2176, a,2 = 0.2082, and 03 = —0.03087. The ratio o" ap (r v )/«v reaches a minimum of 0.65 at c = 4, 
with values of 0.69 at c = 1.8 and 10. Had we assumed isotropy instead, a a . p {r v ) jv v would have been 3% lower (Mauduit & Mamon 2007). 
In this first pass, we assume c = 4, while in subsequent passes we use the relation obtained in ACDM haloes by Maccio et al. (2008): 
c= 6.76(/iM/lO 12 M )-°- 098 . 



6 We used an output of the simulation with 1 particle out of 55 chosen at random. 

7 By T. Beers, see http://www.pa.msu.edu/ ftp/pub/beers/posts/rostat/ 
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iv) We deduce an estimated virial radius r° Bt = ^A/2H with overdensity A = 200. We use H = 100 km s 1 Mpc 1 to 

conform with the positional units of the simulation (and which we had thus assumed for the Hubble flow). 

v) Next, we compute LOS velocity dispersion predicted for the NFW model with anisotropy profile of equation (37) at the projected 
radius of every particle. For this, we use the approximation 



with the coefficients bi given in Table 2 (for the NFW column) of MBM10 and 

GM(r_ 2 )/r_ 2 (In 2- 1/2) c 

vl ln(c+l)-c/(c+l) ^ 

(e.g., eq. [22] of MBM10). 

vi) We filter the particles to have velocities within k times the local LOS velocity dispersion from the global median velocity. We adopt 
K = 2.7, which best preserves the LOS velocity dispersion profile (MBM10). 

vii) We compute the global velocity dispersion of our velocity-filtered sample, this time using the standard unbiased standard deviation. 

viii) We iterate, checking (except after the first pass) that the number of particles has changed or that had changed by more than 0. 1 %, 
by returning to step iii). 

ix) On convergence, we select all particles within 2.7 a z (R) from the median (except for those filtered out by the gapper technique). 

The novelty of this algorithm is that 1) it uses a guess of the (local) LOS velocity dispersion, and 2) it uses a cut at 2.7 a instead of 3 to 
optimally recover a z (R) (MBM10). 




(B7) 
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